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CHAPTER 1 

Purpose of this Work 

Much of what we presently know about "deep space", and the various 
physical phenomena resulting from the interaction of the earth's 
magnetic field with the solar wind, has been obtained from experimental 
evidence gathered by satellites of the EXPLORER series (IMP, for 
Interplanetary monitoring probes) . 

In order to probe, along the same orbital period of a few days, 
the near-earth region, the transition region and free interplanetary 
space, it is convenient to use satellites in geocentric orbits of very 
large eccentricity, typically in the range of eccentricities 
0.9 £ e £ 0.95. Such orbits present a critical "stability" problem. 
Their initially low height of perigee is so perturbed by the graviva- 
tional effects of the sun and the moon that only a judicious choice of 
the launch time can guarantee that the satellite orbit will not exper- 
ience a premature decay in the earth's atmosphere. The determination 
of these times, or "launch window calculation" had very often be 
extremely costly proposition if high accuracy numerical integration 
programs were used. On the other hand, such computation is impera- 
tive when the "target dates" are better known. The present work 
aimed at providing the mission analysts with methods and computing 
tools for studying the stability and evolution of orbits of large eccen- 
tricity. This is the topic tf Chapters 2, 3 and 4. Chapters 2 and 3 



develop an approach for a lower accuracy, but very fast analysis 
technique, whereas Chapter 4 resorts to non-numeric omputation to 
obtain a "symbolic theory", applicable to high eccentricity orbits 
and of average accuracy and computer-time requirements. 

Deep space is also investigated by means of satellites in 
large circular orbits (20 earth's radii, say), which are similarly 
perturbed by the sun and the moon. Although orbital decay is not 
a practical problem here, the development of methods of orbital 
computation, which would be more economical than conventional 
ones, allow for strong perturbations and be singularity-free, appeared 
to be a topic of much relevance. Such is the subjects of Chapter 5. 
Chapter 6 implements the methods of the previous one, and compares 
them to a straight method of variation of parameters, with time as 
the independent variable. 

It is hoped that the approaches and techniques suggested in this 
work will be of help to the mission analyst facing the challenge of 


future space missions. 
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CHAPTER 2 

Method of Approximate Stability Criteria 
2. 1 The Problem 

In the pre-launch phase of the mission analysis of satellites 
in orbits of large eccentricity, there exists a definite need for 
methods of determination of the orbital stability which would combine 
EXTREME COMPUTING SPEED with TOLERABLE accuracy on the results. 

To be more specific, given a launcher of known capabilities, a 
launch site, a spacecraft of (roughly) known mass, an orbit of known 
in-plane geometry (i.e. the initial semi-major axis and eccentricity 
e) , one wishes to "design" an orbit by adjusting, within limits, 

- the orbital inclination on the equator (i a ) 

- the argument of perigee, referred to the equator (m a ) 

- the time (hour; day; year) of launch (which in turn permits 
computation of the longitude of nodes, fi a , and the time of 
passage at perigee, Tp) 

while satisfying, as explained in more detail in the previous chapter, 

- a lifetime constraint 

- other constraints of a technical scientific or operational 

nature (for instance: the angle between the satellite spin 

axis and the earth-sun line, or solar aspect angle, should 
be 90° + 15°, say, at injection) 

Now it should be remembered that the "lifetime constraint" is 
relatively imprecise and can often, to some extent, be relaxed. A 
requirement that the lifetime be "3 years" is not meant to be taken 



2-2 


as a request that the reentry of the satellite Into the earth's atmos- 
phere occur at time t = time of launch +1,095.75 days. Therefore, some 
inaccuracy on the a-priori determination of the lifetime might be toler- 
able when balanced against the speed and economy with which the prediction 
can be made. 

If a method of extremely high economy is indeed obtained, parametric 
studies and the "mass production" of launch windows becomes practical. 
Questions such as this one can be readily answered: "Given this constraint 

on the solar aspect angle, when in the year 1973 should we launch this 
satellite so as to fulfill all constraints? What is the penalty paid in 
lifetime if we make the solar aspect angle condition more stringent? Should 

we try to modify the ascent phase and obtain a different argument of perigee 

11 

etc 

In addition to looking, in this chapter, to a method of appreciably 

reducing the computing time necessary for obtaining a launch window map, 

[ 2— 1 ] 

the accent will be put also, as initially proposed , on automatizing 

the graphical presentation of the output in the form most suitable to the 
users' needs. 


2.2 


Basic Equations 

— - -- y 

Let 0 be the center of the earth, of mass M (Fig, 2,1); r = Os 


the geocentric vector to the satellite; m the mass of the satellite (m is 

infinitesimally small compared to M^, mj) ; m d ,the mass of a disturbing 

body d, assumed to describe a known Keplerian elliptical orbit 

about 0: r,. = r , - Ocl. Now define the vector from satellite to disturbing 

dist d 
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body, p, as (Fig. 2.1) 


-»• -> -*■ 



( 2 . 2 - 1 ) 


The following equations hold, if the vectorial pole 
is the center of mass of the system, 


* ? d 

V +k m d TJ 

d 

R = -k 2 M -^3 


R = -k 2 (M^g- - ) 


for the capital R’s 


( 2 . 2 - 2 ) 

(2.2-3) 

(2.2-4) 


Taking 0 as the origin of vectors, subtracting (2.2-2) from (2.2-3) gives 


r , + k 2 (m, + M) — 3 = 0 


(2.2-5) 


r d 


which describes the elliptical motion of d about 0 (with U<j +M = k 2 (m a + M) . 


Subtracting (2.2-2) from (2.2-4) gives 


r + 




k 2 M -^3 = -k^m(— - ^-j) 
r 3 d • r d 3 P 3 


( 2 . 2 - 6 ) 


Equation (2.2-6) shows that the elliptical motion of s about 0 (with 
p = k 2 M) will be perturbed by the disturbing force due to third-body, m^. 


F 

aist 



(2.2-7) 
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If more than one body is perturbing the orbit of s, or if other perturbing 

forces. F , are acting upon s in a frame centered at 0 and pointing 

other 

towards fixed directions on the celestial sphere, Equation (2.2-6) should 
be complemented to read 


jl. = F +F, F 

r + K m -3 * Aist)1 * £ dist,2 


other 


( 2 . 2 - 8 ) 


in which the F,. . . (i = 1,2,...) have the form given by Eq. (2.2-7) 

GlSt , 1 


2.3 Lidov* s Theory 
2.3.1 System of equations 

In 1961, M.L. Lidov made an important contribution to the problem of 

determining the evolution of satellite orbits under the gravitational 

f2-2 1 

perturbations of external bodies 1 “ . This approach can be summarized 

as follows. Let a, e, i, to, ft be five osculating elements of the satellite 
orbit; p = a(l - e 2 ) is the osculating parameter: v is the true anomaly. 
Angles such as i, to, ft are referred to a plane (such as the equator) 
invariant in inertial space and passing through 0. Along the perturbed 
orbit, the following relations hold 


,dv . dto , 

+ dF + cos 


. dti* f vl/2 
1 -Jj) ■ (up) 


(2.3-1) 


or if 


v = [1 + — F cos 
1 Vie r 


v - — (1 + -)F sin v]" 1 
pe p t 


(2.3-2) 
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r, s, n subscripts indicate components along the radius, the positive 
transverse and the normal direction, r 2 ^ and fiw are related by 


2 dV 
r dt 



(2.3-4) 


A thorough analysis of the order of magnitude of y has been made in the 
course of this grant and is given in [2-3], 

[ 2 - 2 ] 

The results are slightly different from those put forth by Lidov , who 

argued that, since in the expression (2.3—2) for Y, 


F , F 'v — 3 

r v r d 3 


(at most) 


and if e is close to 1, 


r u d r 3 1 

Y " 1±0 ■>, + 


e cos v 


■] £ 


y d a 3 1 

*±<77 7 ^ 

r d 


(2.3-5) 


In the reasoning, an average of r over M, mean anomaly, must be assumed 

since — / 2ir r dM = a. If Lidov' s estimate, as given in (2.3-5), 

2ir o 


rather than ours, is taken to assess the departure of y from unity, one 

(moon) ; 


y d 1 

would conclude for example, that if e = .95, say, and ^ = 8 i 


then if — = - (up to \ of earth-moon distance), or then 

r a 3 3 k 


, . l 1 

1 ±81 27 


20 = 1 + .0091 at most 
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while our estimate (see Ref. [2-3] ) would be 


Y £ 1 + .008 


at most. 


For an extreme eccentricity of 0.95, Y will never depart from 1 

by more than 0.8%. Consistent with the accuracy we are striving to 

[2-2] „ [2-2, p. 720]. 

obtain (of the order of 1%; Lidov quotes 1-3 per cent J» 

and in agreement with Lidov' s approximation, the factor y is taken to be 

1 in what follows. With this approximation, planetary equations are 

written in terms of the five osculating parameters a, i, e, m, 

4^ = 2- /p7 y (e F sin v + F (1 + e cos v)) 

dt e t t 

de _ _ r_ TT/TT v + -JL- 
dt ae ^ t 2ae dt 


dco 

dt 


1 - i £ 

- - /n/u [-F cos v + F t sin v(l + — ) - 
e r t p 


e — F cot i sin(co + v) ] 
P n 


dft 

r 

dt /] 


PP 


sin i n 


F sin(ti) + v) 


di 


— — = /== F cos(co + v) 
dt /yp n 


with e = 1- e 2 . 

Substituting for dt , from Eq. (2.3-4), in which Y = 1, and introducing 

f 2—2 r> 7221 

p = a(l - e 2 ) , the differential system considered by Lidov 1 * ’ is 
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dE . F _ or 4§ , 2 lis (eF + F ) 


dv y t 


dv 


ye 


de _ 
dv y 


(F r sin v + (1 + ~ )F t cos v + e ^ F t ) = - —• 


r v e da 

L p _±-=L 

pae t 2ae dv 


2 

~ = — (-F cos v + (1 + — )F„ sin v - e — cot i F sin(m + v)) 
dv ye v r p t p n 


= — [-F + F_ - sin v] 

ye x t p 


* dn 

cos i ~r~ 
dv 


dfi r 3 1 „ _ . , , \ 

— - = : r F sm(co + v) 

dv yp sin l n 


= — F cos(w + v) 

dv yp n 


(2.3-6) 


2.3.2 Legendre Polynomial (LP) Expansion 

The developments proceed to expand F, the gravitational disturbing 
force, i.e, for body d, 


/P d 


-> -y 
r d - r 


^ xx, + yy . + zz, 

Q\ f '-*■ u ^ p / I d d d ^ 


d 3 

P r d 


3 ' ^'Ir, - *1 3 

1 a 1 


r d ; 


= V M - 

dp r< j3 
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in series of Legendre Polynomials (LP) of argument S - cos 


-*■ 

1 -1 . 
r rd 


-y -y 

rr. 


Specifically, since ~ < 1, with assured convergence 


[ 2 - 4 ] 


.1 _ rr d^ , 

^perturbation ^d p r , 3 


M (r2 + r d ' 2rr d° 


-1/2 r V. 


- ^ [<1 + < f -> 2 - 2< f -) O ’' 72 - f - t ] (2-3-7) 

r, r , r, i. 


. , r 

Now, with a = — 

r d 

(1 + a 2 - 2aS)~ >2 = 1 - | (a 2 - 2as) + |(a 4 - 4a 3 C + 4a 2 S 2 ) 

= 1 + a *S 2 + a 2 (- | + |? 2 ) +.... 

- A«V> 


in which the Legendre Polynomials are 


Po(C) 

= i 



= c 


p 2 (c> 

= - i + -£ 2 
2 2 4 

= |(3? 2 - 1) 

p 3 (?) 

= - 3 C + 5 

2 4 t 2 4 

3 = |(5S 3 - 30 

V° 

■ l' 35 ^ - 

30s 2 + 3) 

P 5 (C) 

= f(63S 5 - 

70S 3 + 15S ) 
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Going back to expression (2.3-7) 


^perturbation ' ^ « l + J 2 


The perturbing force, in turn, is 


F = VR 

d perturbation 


■ i a ^ k "V^r 

r k • - r d' r 

+ V<> 


in which 1 = 3 

r r 


V‘> -3c [ V E » 


A - LSO r • r r ->* r « 

»eV> - • h T - * 0 A cun I 

r d r r d d 


Since ^ = V(r) , the second term vanishes and 


/d li _ t r^- ~ (-) + — — (?) + — — (-)] 
tr ‘ V) ? “ VTT 3x V + r, 3x Sl + r , 3x V J 


■* x d 3 ,x. y d 3 t y\ , Z d 3 ,z^ 

+ V^37 ( r } + ^7i7 ( P + F7a7 ( 7 )] 


We have 


, t r-ii- A + -iL ( i )+ -iL (£)] 

+ z r, 3 z V r. 3 z V r. 3z V J 

d d d 


— (-) 
3x V 


1 _ x_ 3r 
r r 2 3x 



z-10 


9_ ( Z) = _ 3- 
8x'‘r'' r 2 3x 


1.(1) = _ 5- 

Sx^ir r 2 3x 


Gathering terms , 


L. h i , x A * y ^v + ^ 

r d r r d 


+ yy d + ZZ d - 
Vr 

r d r 


= 1 - i-l 

d ^ r r 


Substituting in (2.3-9), and letting q = k-1 


f. - — J 1 (f-) q [[(q + i ) p g+1 (0 - ^ 

d r d ^- 1 r d q+1 


; + 1 ]t r + < + 1 fc)i d ] 


Now, from the recurrence formulae. 


-P - (q + 1)P - £P . , 

q H q q+i 


Finally 


4l + p q« fe)1 d ] 


The components of F, along axes (P,Q,R) are now obtained 


F*Xl 5l Q'^vi asf ^2 


^3 
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The director cosines of the unit vector to the satellite, l^in the satellite 
system (P,Q,R), namely £i> £ 2 * £3 are illustrated on Fig. 2.2. 

In the same system, ^ has components (cos v, sin v, 0). Also 

t, = i o cos v + £2 sin v and let 

2 2 1/ £1 ^2 

5 " <5i + li) » cos = ~ ; sin v = 7 - 

Now, if successive orders in the LP, i.e. q = 1, 2, ..., are considered, 
Equation (2.3-4) gives 

< ? d>, * ~2 <7T>'- p i<'> \ 1 d I 

r d d 


<fy 2 - r? <f-) 2 t- p 2(«) \ I d ' 

r d d 


(2.3-10) 


These are the two values of q retained by Lldov. Rewriting these expressions, 
with pj(?) = 1; P^ (?) = 3s; (?) - y ?2 - \ > 


<*d>. - fl <fr> + r d> 

r d d 


<? d ) 2 = ~^Z 2 [-35 t + |(15C 2 - 3)l d ] 

r d d 


which are Lidov's expressions (4), (5)^ 2 2> P’^ 3 ] # The computation of 
the components of (F^)]^ and (F d ) 2 along (r, t, n) are (Fig. 2.2), since 

rj*r ( r j )* r 

4 = 7“ - — rtr - = 5 cos (v-v ) , 
r j • r r _ *r t> 

a a 


(F.) = [-1 + 3£ 2 cos 2 (v-v )] 

d l.r ^2 r d £ 


(F.) = — 2 7“ f~ 3 ^ 2 COS (v-v ) sin (v-v )] 

’ r d d * * 


< ? d\,n 


~~2 7~ t 3 ^3 COS (v-v )] 
r d d 4 


(2.3-11) 
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and similarly for (Fj) 2 projected along (r,t,n), 

(Fj ) 2 = ^2<f-) 2 [- | ?cos(v-v e ) +y~K 3 cos 3 (v-v ? )] 

’ r d d " 

(F ,) . = “■(— M“ cos 2 (v-v )sin(v-v ) + f? sin(v-v )] 

d 2,t r d r d * 


< f d> 2 ,n=^<t >2[ ^ {2e ’ COs2(V 'V • l ? 3 ! 

d o 


( 2 . 3 - 12 ) 


At this stage, Lidov introduces the notations: 


A = ^- = 1 + e cos v 

def r 


(for the satellite) 


and 


A = — = 1 + e cos v j (for the disturbing body) 

d def r a d 


after which the components of the forces are rewritten: 


<Vi, 


= 3 -^5- — (-g 6 + cos 2 v + 2P3 sinv cos v + 02 sin 2 v) 


r ' Pd 2 p d 


& 


y p 

) = -3 -S_ — f sin v cos v (0* - 62) + (sin 2 v - cos 2 v )03]^ 

lft Pd p d 


iA 1 

(F,), = 3 -^2 ~ (& 5 cosv+ sin vV 

d l,n p/ p. ' 


Pd F d 


(F,V r = ¥~ ^2 | <*1 cos v - | a 2 sin v + y l cos 3 v 

d 2,r Pd p d 


• + 3 y. cos 2 v sin v + 3 Yc cos v sin 2 v + y sin 3 v)JL 


2 ' 2 

A 


oil- 1 
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y 2 

(F d ). t = ^*^2 ^T[-Y 3 cos3v “ 2Y 6 )cos 2 v sin 

u ^ D , D, 


V p d 


~(2y - y 9 )cos v sin 2 v + y 


sin 3 v + - a, sin v 


1 

- a-2 cos 



V 2 

(F,)„ = ^ — tOiV cos 2 v + 2y 7 cos v sin v + y 6 sin 2 v 

d 2>n 2 Pd Pd 


in which 


5 V 


(2.3-13) 


a l “ 5l Ai» a 2 “ ?2 A j» 


«3 = ?3 


2 3 
A d* 


2 3 

= 5 2 V 


^3 = e^V 


^4 

y i 

Y4 

Y7 


^2 ? 3 A d* 
3 4 
^l A d’ 


5 l ? 3 A d* 


«lWd’ 


3s 

Y2 

^5 


M3 A d * 

3 4 

S 2 A d ’ 


S 2 ? 3 A d’ 


3,6 

V3 

Y c 


- A 


d * 


Wd* 


«2«lV 
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d 3. ^|p 

Going back to system (2.3-6), in which is taken rather than in 
the first equation, and substituting for the forces, their expressions 
in ( 2 . 3 - 11 ), ( 2 . 3 - 12 ) etc., r.h. sides are obtained which depend on 
the motion of the disturbed body (v) and on the motion of the disturbing 
body (v d ) through 


- powers of r = 



P 

1 + e cos v 


- positive powers of sin v and cos v 

- the greek symbols a., B.» Y, etc. of the form 

1 J K. 

(2.3-u) 

1 2 3- k 

in which Ci , C 3 are the director cosines of the unit vector 

4 - 4 4 . 

to the perturbing body in the satellite system (P , Q, R) . 


It is taken for granted that, consistent with the approximation in a 
first-order theory, the satellite orbital elements a, e, i, etc., are 
taken as constants in the r.h. side of system (2.3-6) and that the 
changes in these elements are computed separately, for each disturbing 
orbit, over one orbit of the satellite and then linearly superposed. 

To summarize: suppose that we have a suitable representation 

3 b c ^ 

for Ci, ? 2 > 5 3 , 4^, given the true anomaly v of the satellite. Then, 
for any element z for which an equation in system (2.3-6) is written 
(the set of elements z is denoted z) 
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After integration of (with respect to v) from 0 to 2ir , where subscript 
I is the order of the LP expansion, i.e. I corresponds to forces 


* 0 [ (— ) 

ra 


1+1 r 2+1 

], II to forces M3[( — ) ] etc. 


we obtain 


ic 2 

(Az) t = change in element z due to forces of M)[(— ) ] 

I r A 


and, in total. 


(Az) = 2 (Az). 

j = I, II, III. .. . 3 


(2.3-15) 


(2,3-15) expresses the change in any orbital element due to the various 
orders in the LP expansion, i.e. ordering these in columns 


LP 


(Az). 


(Az) 


(Az) 


II 


III 


(2.3-16) 


2.3.3 Taylor Series expansion 

ft 

Lidov’s theory then further proceeds to expand any of the above re- 
ferred greek symbols a^, , 8^ 5 y^ as 


act. 

a i = (a i ) ref + ( dt^ ref At + ( ' 


dV 

^ref 2! 


At' 


(2.3-17) 


in a Taylor's series about a given point in time, t along the satellite 

revolution. In Equation (2.3-17), At is the difference t - t re f» measuring 

the time elapsed from t ,. From (2.3-17) and (2.3-14), it follows that 

ret 
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for any (Az) I# (Az) ri , ...» any of the rows of (2.3-16) can be further 
divided as 




— * Taylor Series 



< 4z) l,l 

(Az) 1,2 



<4z) 2 _i 

(4z) 2,2 


LP 

> Wz) 3,l 

<Az) 3,2 

(2.3-18) 


thereby providing 2 directions of expansion: the first one (i, first 

vi+K 

subscript) corresponds to the LP expansion of the forces (term ^ ) ) > 

d 

and the second one ( j , second subscript) to the number of terms, or the 
order increased by one, retained in the Taylor Series expansion of the 
quantities , 8^ , .... 

To appreciate what is involved in the integration of (Az)„, we 
shall look specifically at the cases: 

a ) j = i 2 , 3 i = 1 (FIRST COLUMN IN TABLEAU (2.3-18)) 

This specifically amounts to assuming that the disturbing body 

is f ixed in space during any one revolution of the satellite . (The 

disturbing body is "frozen" at an average position) One should 

expect this approximation to be the better the smaller the ratio 

of — (or — ), since the a’s, S's, etc. would indeed be sensibly 
r d r d 

constant if the period of the satellite was infinitesimally small 
compared to the period of the disturbing body in its apparent geo- 


centric motion. 
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The terms to be Integrated will involve 
bed 


k h h A d 

a b c 

appearing as constants: (£ 1 ). (S?),. * (5*)*. ( 1+e d cos v d^ref » 

^ref c re£ c ref 


and At does not appear , 


b) 1 = 1, 2, 3, ... , .1 = 1, 2. (FIRST AND SECOND COLUMNS IN TABLEAU 

(2.3-18)) 

Here it is of course assumed that the oc f s, 3 1 s, etc. are suffi- 
ciently well described by a straight line tangent to the corresponding 
curve a = a(t) at t = t re f, (Fig, 2,3). This approximation should 
hold well if the angular motion of the disturbing body is "slow 11 
compared to that of the disturbed body. For our purpose, such an 
approximation is amply sufficient for the sun's contribution to the 
perturbations " . The additional terms to be inte- 

grated (compared to a)) will involve expressions originating from 
abed 

and reading like 

'v At x[constants computed at t = t^^] 

Similarly, the third column in tableau (2,3-18) accounts for terms 
in (At) 2 etc. .. 

The assessment of the order of magnitudes of each contribution 
(Az)_ has been done in detail in [2-3] 


2,3,4 Results of Lidov's theory: short-range and long-range 

As an example, consider the "11" theory. Namely, only forces of 
r 2 

0(— ) are retained in the LP expansion, and furthermore, the perturbing 
r , 
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bodies are fixed at the position they assume at t — ! re f* Taking for 

di 

instance the equation for , 

•li = — F cos(oH-v) 
dv yp n 

Replacing (F n ) by (F n ) ^ as given in (2.3-11), 


y d » r 


= 3 — (— ) “ cos(w + v)[?. cos v + 5, sin v] 

dv y p 3 i z 

3 

- 3 „ A d 

= 3 — (— ) e 3 £~— 4 cos(w + v)[C 1 cos v + sin v] 

V Pa 3 A 1 


It is apparent that the evaluation of the following definite integrals is 
required 


. 9 tt . ^ cos^v . ,2*rr,^_ *_ v sin v cos v 

{ ( V^3> — dV * { (A d^ 7 — dv ■ 


.2tt _ v cos v sin v 

/ 3 dv > 

u A 


/ 3 sin 2 v 

o (Aa^s^r dv • 

A 


One suspects that as the indices i and j are increased beyond 1, the 

volume and complexity of the calculation might become prohibitive. Hence, 

resort was made, in the present project, to non-numeric manipulation on 

the computer for the development of a modified, extended Lidov's theory. 

This is treated in Chapter 3. 

r 0—2 1 

In his paper , Lidov gives for the five elements a, e, i, w, ft, 

the results of the "11", "12" and "21" theories in tableau (2,3-18). 
Limiting ourselves to the "main" contribution "11", we reproduce Lidov's 
results (SHORT-RANGE PERTURBATIONS): 
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A^a = 0 


A n e = -!5it ~ (™) ee 1/2 g 3 


A n Q - 15tt (J-) 3 “n^ [(1 - \ e)3 5 sin u 

M p d e 2 sin i 


1 

+ - e( 5 4 cos to] 


A n i = 15tt ^ (|-) 3 [ (1 - | e)S 5 cos to - ^ e3 4 sin to] 


d e 


p d a 3 1/2 

A n to = 3ir — (— ) e [43x - 02 “ M ~ A ll n cos 1 ( 2 ‘ 


y P, 


In the last part of his paper, Lidov investigates the secular 
changes in the elements of the satellite orbit by integrating the orbit- 
to-orbit changes, as given by Equation (2.3-19), over the period of the 
disturbing body . He thereby obtains the following expressions for the 
secular changes $z in the orbital elements z, per satellite orbital 
period 6t (LONG-RANGE PERTURBATIONS) 

S3 L 

6 ns. = 0 

A ll e = ^ A^ee sin 2 i sin 2w 


3-19) 
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6 1 


- - I k d ^ITT t(1 - E) sln2 “ + l E) 


(2.3-20) 


5 11 U 


2 

2 A d “T/2 [(cos i_e ) sln2u + 5 


In the above equations, the plane of reference for measuring the 
angles is the orbital plane of the perturbing body d , and for disturbing 
body d, is defined here as 


. V d a . 3 3/2 

d U p d d 


(2.3-21) 


Finally, using Equation (2.3-20), Lidov is able to classify the 
long-range behavior of the perturbed orbits in terras of the two integrals 
which, besides the trivial one: a = constant, could be determined, namely 

ci = (1 - e 2 ) cos 2 i = ecos 2 i 

c -2 = (1 - e) (| - sin 2 i sin 2 to) (2.3-22) 

These integrals, which apply to the system of differential equations 
describing the secular change of the orbital elements due to one perturbing 
body, in the absence of oblateness, will be used in the approximate 
stability criteria method which follows. 
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2.4 The. Approximate Stability Criteria Method 

2.4.1 Introduction 

The approximate stability criteria method was developed under this 
grant by Renard^ 2-5 ' 2-6 ^. Its goal is to provide a fast, economical 
(if less accurate) way of determining the stability of an orbit of large 
eccentricity and, in final analysis, the quick generation of launch window 
maps called for in a mission analysis. The method has since been used 
with success to study the launch windows of several satellites of the IMP 
(Explorer Series) t 2-7to 2-9 and is operational at NASA Goddard Space 
Flight Center. 

2.4.2 Some definitions: 

_ Orbit of large eccentricity : this is defined here as a geocentric orbit 

having an eccentricity in the approximate range 0.9 £ e £ 0.95, or equi- 
valently a geocentric distance to apogee R^ £ 20 to 40 (earth's radius), 
if an initially low perigee, close to the earth's surface is assumed. 

- Stable orbit : rather than being called (as it maybe should) "successful" 

an orbit is called stable if the height of perigee h^, remains during the 

* 

whole spacecraft lifetime, L, larger than some critical value h^, equal to 

& 

or slightly lower than the initial value h^. h^ corresponds to an assumed 
height of perigee leading to orbital decay in the atmosphere. 

— Launch window, launch window map : A launch window is the set of points 

DL, HL (day of launch, hour of launch) for which stability is realized, 
and for which a number of technical, scientific or other constraints are 
met^ 2 ” 10 t0 2 “ 16 ^. The boundary of the launch window defines the so-called 
launch window "map" (Fig. 2.3). 
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2.4.3 Present method and criteria 

2. 4. 3.1 Evolution of the orbital elements 

Fig. 2.4 (a. to f.) illustrates the evolution with time of some 
characteristic quantities of high eccentricity, stable orbits: the 

altitude of apogee, the altitude of perigee, the inclination of the 
equator, the longitude of nodes, the eccentricity and the argument of 
perigee (Ref. [2-5]). It is noted that for such a stable orbit, and 
a dense satellite, due to the rapid increase in perigee height, the 
effects of Earth’s oblateness air drag are very limited and affect 
stability rather indirectly. Thus, in first approximation, they will 
be neglected in the analysis. Adjustments to the lifetime estimation 
might have to be made, however, in those special cases where the 
effect of oblateness plays a more significant role, as is mentioned in 

Chapter 3. 

2. 4. 3. 2 Motivation 

A purely numerical determination, on the computer, of the launch 
window map for a satellite having a required lifetime of at least one 
year could require an average of 50 to 100 hr. of IBM 7090 per year 
of possible launch dates. Addressing herself to this problem, 

M. Moe^ 2 developed simplified equations which were later solved 

on the analog computer^ 19 >2 19] a j- a considerable gain in computa- 
tional speed and with good agreement between the predicted and exact 

, [ 2 — 10 ] 
window contours 

It remained tempting, however, to try and define the launch window on 
the basis of approximate stability criteria which if fulfilled at launch. 
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would very likely guarantee the whole lifetime . This would result in 

an economical and fast method on digital computers, more universally 

available, and allow for better ephemerides of the Moon and the Sun, 

the orbits of which were taken to be circular in the above mentioned 

f 2-2 0 1 

papers. Leroy and Pace 1 had mentioned Lidov's theory (Ref. [2-2]), 

as a possible way to somewhat restrict the domain to be investigated 
numerically. With the goal of establishing stability criteria, we 
found very encouraging that many of the orbital features just described 
were qualitatively predicted by Lidov's analysis. Of course, the domain 
of validity of Lidov's results was presumably restricted to lower 
eccentricities than those retained here. For example in Ref. ([2-2]), 
Lidov was aiming for an accuracy of 1 to 3% with geocentric orbits of 
semi-major axis of the order of 30 to 40 x 10^ km. These figures were 
perhaps too conservative, since analog integration of the similar 
M. Moe's equations had given good predictions of the launch windows, 
up to eccentricities of the order of 0.95. 

2.4.3. 3 Setting up the criteria 

For stability, we require that r p , radius at perigee, not decrease 
with time, over the satellite lifetime L. Let a, e. be the semi-major 
axis and eccentricity at perigee, respectively, and 6z the change in 
quantity z from one perigee to the subsequent one: 

6r p = 6 [a (1 - e)] = 6a- (1 - e) - a6e (2.4-1) 

According to Lidov's "11" theory (and this is also true for the "kk" 
theory, k > 1, see Chapter 3) 


6a = 0 
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Thus , 

6r = - a<5e 
P 

Thus, it is required that the orbital eccentricity e have a decreasing 
trend with time. The constancy of a, as obtained from computer results, 
is illustrated in Fig. 2-5, for the orbit of Fig. 2-4. 

The principle of the present method is to simultaneously require 
that stability be realized: 

1) In the long-term (subscript LR, for "long-range"), having 

characteristic time T^(moon) or Tg(sun): CRITERION 1. 

2) In the short-term (subscript SR, for "short-range"), having 

characteristic time t , i.e. a few days: CRITERION 2. 

Set L 

3) In the intermediate-term, so the "waviness" of the curve of 

height of perigee vs. time about its trendline is limited 
(characteristic time x^/ 2 or x Q / 2): CRITERIA 3, 4, 5. 

4) In the very-long term (characteristic time T vUR ’ as ^ et 

unknown): CRITERION 6. 

These various stability criteria are now studied one by one. 

LONG-TERM STABILITY 

In the long-range, stability should exist for the secular effect of 
Sun and Moon, i.e. on the average over a period of the perturbing body. 
As obtained by Lidov (Ref. [2-2]) for one perturbing body "d", and re- 
called in Equations (2.3-20) above: 

6e,,= y A. e e sin 2 i, sin 2(i>, 

1 1 4 d d d 


A 


d 





(2.4-2) 
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Note that i,» w, are the inclination of the satellite orbit, and the 
d a 

argument of perigee of the satellite, referred to the orbital plane 
of perturbing body "d" . ' 

For the sun and the moon acting simultaneously, we can therefore 
state the long-term stability criterion (CRITERION 1) 

6e^ R = ^e e ^2(A q sin i^ sin 2 u)q + sin i^ sin 2(0^) - 0 (2.4-3) 


This will define a long-range stability region 
plotted in terms of launch hour vs. launch day. 

It should be noted that the ratio of the 
and Sun, respectively, is 

« 2.18 

A a "q p M e ® 


, which can easily be 


amplitudes A^, for Moon 


and is obviously independent of A. As an example, for an orbit having 
the following characteristics: 


h s = 203,632 km 
A 

h - 192.6 km 
P 

e =0.93932 

t = 4.1047 
sat 

we obtain 


6h 


p ,LR 


6h 


p,LR 


max 

max 


due to the Sun = 51.9 km/revolution 
due to the Moon = 113.2 km/revolution 
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Computer studies of high' eccentricity orbits have shown that, starting 
with a relatively low inclination 0 < i < ^ on the ecliptic, stable 
orbits were accompanied by a significant increase in inclination on the 
plane of the disturbing bodies. It should also be noted here that if 
only one predominant perturbing body is considered, or at high inclina- 
tions on the ecliptic, in which case ig % i^,the qualitative use of 
Tisserand’s criterion made in [2-11] to account for the eccentricity vs. 
inclination relationship just appears as a quantitative consequence of 
Lidov’s secular theory to order "11". Indeed, 

= - | A(1 - e) sin 1 cos 1 sin 2(0 (2.4-4) 

Dividing (2.4-4) by (2.4-3), and after some manipulation 

1/v 

6 (cos i) _ _ 6e 

cos i "" 1 /2 

e 

or 

(1 - e 2 ) 1 ^ 2 cos i = constant (2.4-5) 

which is one of Lidov's "secular" integrals. Obviously, for stable 
orbits of j < i < tt , the inclination will decrease as time increases. 

SHORT-TERM STABILITY 

As recalled in Equation (2.3-19), Lidov’s "11" theory gives for the 
short-term (subscript SR, for "short-range") change of e due to dis- 
turbing body "d", and per satellite revolution: 
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S^e = - ee 


>/ 2 A 


3. 


e d 


3 /2 P 3,d 


(2.4-6) 


3 


in which 6 = C 1>d C 2>d (/) rer Fi 8‘ 2 * 2 shows the geometrical 

significance of B . If the projection of r d on the satellite orbital 
plane is in the FIRST or THIRD quadrants, then 8 is > 0, and from 
Equation (2.4-6), the orbit is stable in the short-range; if the pro- 
jection is the SECOND or THIRD quadrants, the orbit is unstable in the 
short-range. 

Now, for the two disturbing bodies (Sun and Moon), we require 
short-term stability by stating CRITERION 2: 


Se SR ' 


~ ee 


/ 2[( A 


5/ 2 > e 3,M +( ; 

£ M E 0 


0 


3 / 2 ^ 3 3 , 0 ] 


(2.4-7) 


as shown in Fig. 2.7. An alternative form of CRITERION 2 is 

<^ 3/2A «3,M + ¥,,. ; ° 


(2.4-8) 


INTERMEDIATE-TERM STABILITY 

Even if the short-term evolution of the eccentricity is favorable 
initially, 3 3 M will change sign as 1^, unit vector to the moon, 
rotates in inertial space. Therefore, the eccentricity will oscillate, 
over the lunar month, about an intermediate trendline , which corresponds 
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to the long-range effect of the Moon and the short-range of the Sun, 
averaged over x^. The slope of this line should be < 0 (Fig. 2-8), 
which is expressed in CRITERION 3: 


6e lNT =<$ e sR.0 > T + (Se) LR,M <0 
M 


(2.4-9) 


If the latter condition is fulfilled, it remains to require that 

the waviness of the eccentricity vs. time curve (or alternatively, h^ 

vs. t) not be so pronounced that e increases again towards its initial 

value, at its next maximum (or h decreases again towards its initial 

P 

value, at its next minimum), j designating the index of the perigee 
passage (the initial perigee has index 1) corresponding to that mini- 
mum in 11^, the above requirement is approximately expressed by CRITERION 
4 (or 2+3 strong) . 


6e 


j : x , A M „(k) . A a Jk) 


j = e j “ e * “ k=0 [ F l/ 2 3 3,M + 3/2 e 3,0 ] 


2.4-10) 


'M "0 

It is apparent that criterion 4, which limits the tolerable lunar 

modulation, encompasses criteria 2 and 3, but these are taken in this 

order because they are more readily checked than 4. Criterion 4 is 

then disregarded if 2 or 3 leads to a failure. 

Now, the same reasoning is repeated for the solar modulation of the 

eccentricity vs. time curve (Fig. 2-5). The criterion corresponding 

here to 3 is 6 e <0 (criterion 1), and the one corresponding to 

criterion 4 is now developed* The upper limit on 5 e on account of the 

LR 


solar modulation about its trendline may be simply approximated, near 



2-29 


2 

the limit of stability,- by 0 >T “ 5 e LR l ~ ) (amplitude of 


sat 


wave 


>' * 972^ (tlme t0 C ° ntaCt) ’ l^f ' Se <SR,8>T M - Se VS? 


CRITERION 5 (or strong) is then stated as 


(a) 


6e LR < 9/3 ^ 5 e <SR,0>T M 6 e LR ) lf 6 e <SR,0>x M < 


<5e 


LR 


(2.4-11) 


(b) taken to be satisfied if ^ e <SR q> > < ^ e x J R (2.4-12) 

t m 

All quantities in (2.4-11) or (2.4-12) have been determined previously. 

VERY-LONG-TERM STABILITY (FOR ORBITS WHICH ARE NOT QUASI-NORMAL TO THE 
ECLIPTIC) 

It remains to ensure that the very-long range effect of the 
motion of the apsidal line of the satellite orbit with respect to the 
orbital planes of the perturbing bodies will not cause the eccentricity 
to reach its minimum value before half the expected lifetime has elapsed 
(Fig. 2.9). 

Computational results, for example those of Fig. 2.9 for IMP-C 

under the solar and lunar influences, suggest that a max " a might be 

approximated in the region of interest, and when the inclination of the 

orbit on the ecliptic or the moon's orbital plane is not near 90 ° 

(more is said about this in the next section) , by a half-sine wave with 

unknown "very-long-period" x ±n (Q> As is shown in F ig. 2 - 10 , 

V LK 

e - e = (e -e.) sinir ~ — (2.4-13) 

max max min T VLR 

Assume that there exist one predominant disturbing body. The non-trivial 
integrals in Lidov's "11", secular theory are, as given before, 



2-30 


C£ — (1 - e 2 )cos 2 i = ecos 2 i 

C 2 « (1 - e ) (~ - sin 2 i sin 2 co) (2,4-14) 

From these, the extremal values of e = 1 - e 2 can be determined from: 


a) 


if C 2 > 0 


e 


max 



^2 


e . 
mm 


|[1 +| (c x + c 2 ) - { (1 + | (ci + C 2 )) 


2 


- Fc.) l/2 1 


(2.4-15) 


b) if c 2 < 0 

e , e . are roots of the quadratic equation 
max mm 

e 2 - [1 + - (c ^ + c 2 )] e + ^ cj = 0 (2.4 t16) 

Now, in order to be able to compute c^, c 2 , for a given initial orbit, 
we should have a unique plane of reference (orbital plane of the dis- 
turbing body) , with respect to which angles i and to of the satellite 
orbit are measured. For moderate inclinations i on the moon's orbital 
plane, in view of the small value of the inclination of the moon's 
orbital plane of the ecliptic (i = 5.145°), and of the dominant effect 
of the moon, it can be assumed that the sun approximately describes an 
orbit coplanar with the moon's orbit (this approximation would break 

down if i £ y , and for a perigee located "between" the moon's and the 

[ 2 - 6 ] 

sun's orbital planes ). 
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In this approximation, ci and C£ are computed with values of 
± 9 a) referred to the moon's orbital plane . This is the approximation 
adopted in an early version of the SABAC program, and illustrated in 
the example of Section 2.5. 

For more accuracy, we later decided to use values of and C 2 
which would account for the unequal magnitudes of the amplitudes of the 
perturbations ("A^") , for sun and moon, by retaining as values of Cj 
and C2 the weighted averages 


c . 
3 


A M + A 0 


(c j Hi only 




0 only 


0 - 1 , 2 ) 


(2.4-16) 


in which (c^) and ( c j) 0 only are computed from (2.4-14) for angles 

(i w„) and (i , to ) , respectively. This procedure is the one em- 
MM 0 U 

bodied in program SABAC (Version A) . 

So far, t is an unknown quantity. If the orbit has been found 
VLR 

to be stable in the intermediate range (criterion 3 satisfied), one can 
use Equation (2.4-13), to define an angle (0 < ~ 2 } ’ w ^ ere t ^ ie 

eccentricity is e 0 


sin = 


c max 


- e c 


e - e . 
max mm 


0 < e < ; 


6e 


The slope at £ = K is estimated to have value 


riSW 1 ) for 

sat 


t of the order of one year: 
VLR 


<5e 


INT 1 


sat 


<5e f \ e 

— - = ( e - e . ) »cos 
6t max min 


VLR 
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Normalizing to t , we obtain t vlr , 


t vlr 


7TT ^ (e 
sat max 


e . ) cos 
mxn 


Co 



(2.4-17) 


* > £ 0 

The lifetime condition is satisfied if T = t (1 - 2 — ) is 

V-LK 

larger than L, required lifetime. Thus, CRITERION 6 is stated as 
follows : 
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VERY-LONG RANGE STABILITY (for orbits quasi-normal to the ecliptic) 

When the inclination of the orbit on the moon's orbit (or equi- 
valently, on the ecliptic) is in the neighborhood of 90°, the lifetime 
criterion (6) is modified as follows, for those orbits whole perigee 
motion, during a significant fraction of the lifetime, occurs between 
the orbital planes of the two perturbing bodies (Fig. 2.11), for 
instance under the ecliptic and above the moon’s orbital plane. 

Since u-w is then very small, it is no longer valid, even though 
^ i that ^ This explains why the above described criterion 

leads to predicted lifetimes which are systematically in excess of 
actual values. 

It was found that typically, for an orbit of the IMP-G type 

(5 0 

(e c ~ 0*946; % 90°), — was very small and e would vary over a 

one-year lifetime between 0.946 and 0.934. Therefore, assuming that 
the sun describes an orbit in the moon’s orbital plane amounts to 
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neglecting the unfavorable effect of the sun on the (already small) 
change in eccentricity due to the moon. 

The alternative approach taken in the present case (and imple- 
mented in digital program SABAC version B) is as follows. A plane to 
is defined, which is obtained by rotating, about the nodal line of the 
moon in the ecliptic, the plane of the ecliptic towards the moon’s 

a m 

orbital plane by an amount i^g , where i^g is the mutual 

inclination of these two planes (5.145°). For the theoretical 

justification^ 2 ” 6 ^, see the treatment contained in Chapter 3. 

Now, let to, be the argument of perigee of the satellite orbit, 

referred to plane to. The eccentricity will reach a minimum when the 

perigee will be exactly contained in plane to, after a time equal to the 

half-lifetime, 1/2. For confirmation, the reader should refer to 
* 

Chapter 3. For small sin to, in the case of a southwards injection, 
we write CRITERION 6 (Version B) : 


_ tt - to, > !L 
2 ,6to, , 2 


(2.4-18) 


The time-rate of change of to, is computed from + ^6 T ^LR,S^ + ^ ’ 

In this formula, from Lidov's formula in the long-range 


(— ) 

V ’6x LR,d d 


= A d [(cos 2 i d - e)sin 2 to d + | e]~^- 

2e 


(2.4-19) 


Factor (1 + w) results from averaging the expression between brackets 
in C|“) SR > due to the sun » i<e * > with e Q ^ 0 , 


.6to 


A 0 2 


.0(0, -o _2. £ '2rA £ _ r _ 1] 

'6t^ SR,Q = 5 1 n,0 ^2,0 J 


(2.4-20) 
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The average is taken over half a solar modulation (t^/ 4), to obtain the 

2 

linear trend of uL vs. time, and £[ can generally be neglected in 

* 

Equation (2.4-20). More details and examples are given in Chapter 3. 


2.5 Examples of Application of the Method of Approximate Stability Criteria 
From 1968 to the date of this writing, the approximate criteria approach 
has been used as a fast, economical tool to generate the launch windows 
of satellites in orbits of large eccentricity. The present section will 
deal in some detail with examples of application on satellites: IMP-B, 

IMP-G, IMP- I, IMP-K and K' (mother-daughter system). 

2.5.1 A check of the method: IMP-B launch window 


In order to check the effectiveness of the above method, it was 
decided to try and recover the launch window map for IMP-B, which had 
been well documented^ 2 1 1 ^ . 

This launch window map had been established during the preliminary 
studies on IMP-B and -C. The orbit had the following initial data 


(at injection): 


= 109,952.5 Nmi 
» 104 Nmi 

= 108,290.5 km 
= 0.93932277 
= 32.912693 deg ” 
= 133.659044 deg _ 


referred to a (earth's equator) 


Days studied: April 11 to June 15, 1965 

Hours studied: 8.00 to 18.00 hr. U.T., time of injection, at perigee 
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Failure to meet any of the criteria, 1 to 6, led to the rejection 
of corresponding launch hour on the day considered. The results are 
compared to those obtained by NASA's numerical integration program ITEM, 
based on Encke's method (Fig. 2.12). It is seen that the topologically 
complicated features of the contour separating the "stable" and "unstable" 
regions are well recovered. The error in predicting "peaks" or "valleys" 
in the contour is at most of the order of a few tenths of an hour, and 
much less on the average. The largest discrepancies are recorded at ITEM 
"marginal" points, i.e. for orbits, otherwise stable, having a height of 
perigee between 90 and 100 N.mi. for one orbit only, which is of little 
consequence in practice. This accuracy appears sufficient for the 
purposes of mission analysis. 

It should also be emphasized again that the lifetime condition 
is far from being the only one to be considered. Constraints of a 
technological or scientific nature will further reduce the stable region 
into a much smaller one, acceptable for the mission. In this reduced zone 
of the maps, a final, accurate sutdy is then made, using elaborate and 
expensive digital integration methods. 

2.5.2 General comments on the economy of the method 

It is obviously impossible to accurately pinpoint the savings 
factor obtained by using one method compared to another, particularly 
in a time-sharing environment. However, good estimates of the orders 
of magnitudes can be given, and to the maximum extent possible, the 
conditions in which comparisons are made will be clearly stated. 
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As an implementation of the above method, a program called SABAC 
(for Stability Analysis By Approximate Criteria) was written in FORTRAN 4 
to check, on the basis of the above criteria, whether a point on the 
launch map is stable. 

On the UNIVAC 1108 of Carnegie-Mellon University and in OS, it 
took no more than 0.02 sec per calculation point. This figure should be 
compared (with, as we said above the "order-of-magnitude" viewpoint) to 
7 to 10 min for conventional numerical integration programs run on the 
same machine, integrating over a one— year lifetime. Hence, an economy 
factor of the order of lo\ with the same amount of information obtained 
in the determination of the overall launch opportunities. 

It is worth mentioning that SABAC includes an analytically defined 
ephemeris of the Moon, giving the distance with an error smaller than 
500 km at maximum. The Sun's ephemeris is read in. 

As a last comment, it should be repeated that the method is obviously 
no substitute for the detailed study, by a numerical integration on a 
digital computer, of a particular set of launch days and hours, and the 
corresponding history of the orbital elements over the whole lifetime. 

But this may now be done only in those finally selected "target" regions 
of the map, where all conditions of constraint are met. 

Program SABAC comprises 2 versions, which differ by the method used 
to estimate the lifetime. Version A is suitable for orbits which are 
not quasi-normal to the ecliptic, i.e. it should not be applied to IMP-G. 
Version B is suitable for orbits nearly normal to the. ecliptic , such as 
the orbit of IMP-G. Both versions are thoroughly documented as part Dl 
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of the volume "Documentation of programs and subroutines” appended to this 
report. A slightly different version, of the same programs, also exist 
at NASA, the difference being that a subroutine, SUNEPH, computes the 
Sun’s ephemeris rather than entering the Sun’s coordinates as data, 

2,5.3 IMP- I launch window 

IMP-I (IMP6) was a 636 lb. spin-stabilized spacecraft, with its 

[ 2 - 9 ] 

spin axis nominally perpendicular to the ecliptic plane . It 

carried a payload of 12 scientific experiments and one engineering ex- 
periment. It was launched on March 13, 1971 at 11.15 EST, and inserted 
into an orbit having the following characteristics 

- Orbital period : 4.13 days 

- Perigee : 243 km (initial) 

- Apogee : 206,258 km (initial) 

- Lifetime in orbit: 3.6 years 

- Inclination : 28.69 deg. 

A preliminary orbit was given to us by GSFC, in late 1969. It 
had the following parameters Revised, 02/70 

- Height of Apogee : 217572.19 km 216676.62 km 

- Height of Perigee : 277.7998 km 240.24 km 

- Inclination (equatorial): 28.90053° 28.2996 

- Argument of Perigee: -53.1456° -66.2037 

- Longitude of Perigee: 115.91055 (East) 112.67 E 

Latitude of Perigee : -22.75° (South) -25.838 S 

- Lifetime : 3 years 3 years 

The launch windows of Fig. 2.13-2.14 were obtained. 
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It is of importance to note here that a NASA GSFC’s request, a 
series of cross-checks were made between NASA’s digital integration 
program (ITETCd) , SABAC and C-MU’s program EOLA. Excellent agreement 
was found, provided the allowed drop in perigee (73 km) was kept in 
mind when inputting the data of SABAC. The results are summarized in 
the following table. 
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Cross-check of ITEM, SABAC-1A and EOLA 



Day* 

TIME* 

ITEM 

SABACl-A 

E0LA 

1) 

319 

1 

Success** 

Success 




2 

835 days 

777 days 




18.5 

Success 

Success 

Success 

2) 

320 

19 

Success 

Success 




20 

Success 

Success 




21 

Success 

Success 


3) 

324 

20 

Success 

Success 


4) 

325 

19 

Failure on 4-th Or- 
bit; perigee ht. = 
187 km. 

Ripple fail. 
(Criterion 4) 

Perigee history; 
240, 267, 168, 
327, 189, 142 km. 

5) 

326 

19 

Failure on 4-th Or- 
bit; perigee ht. = 
134 km. 

Ripple fail. 
(Criterion 4) 

Perigee history; 
240, 228, 182, 
288, 134 km. 

6) 

328 

19 

Failure on 1-st Or- 
bit ; perigee ht . = 
149 km. 

Short Range 
fail (Crit. 2) 

Perigee history; 
240, 149 km. 



20 

Perigee ht. drops 
to 190 km. in 4 or- 
bits (drop of 29 km. 
in 1-st orbit). 

Short Range 
fail (Crit. 2) 

Drop of 27 km. 
in 1-st orbit; 
drops to 200 km. 
in 4 orbits. 



21 

Success 

Success 

Success 

7) 

330 

19 

Failures on 1-st Or- 
bit; perigee ht. = 
153 km. 

Short Range 
fail (Crit. 2) 

Perigee ht. drops 
to 148 km. in 1-st 
orbit. 



21 

Perigee ht. drops to 
199 km. in 1-st 
orbit. 

Short Range 
fail (Crit. 2) 

Perigee ht . drops 
to 202 km. in 1-st 
orbit . 

8) 

331 

21 

Perigee ht. drops to 
161 km. in 1-st 
orbit. 

Short Range 
fail (Crit. 2) 

Perigee ht. drops 
to 163 km. in 1-st 
orbit . 


* 

Day: 
Time : 
**Success : 


Day number, 1970 (Reference - Jan. 1=0). 

U.T. in hrs. at injection (assumed to be at perigee). 
Based on a 3-year lifetime. 



DAY* 

TIME* 

ITEM 

S ABACI- A 

E0LA 

332 

22 

Perigee ht. drops to 
186 km. in 1-st 
orbit . 

Short Range 
fail (Grit. 2) 

Perigee ht. drops 
188 ton. in 1-st 
orbit. 

334 

24 

Perigee ht. drops to 
202 km. in 1-st orbit . 

Success 



SABAC1-A results reported here permitted a perigee drop of 73 ton. (DPLIM = 73 ton). 
If no perigee drop were permitted, the failure points would remain so, 
while a few of the "success" points would turn into failure. 


E0LA is a digital integration program based on a variation of parameters 
method with the true anomaly as the independent variable. Earth's oblate- 
ness was included in these runs but the atmosphere was not. 
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B. Kaufman and D.P. Muhonen, at NASA GSFC, carried out a detailed 

r 2-9 "1 

analysis of the IMP-I launch window 1 , subject not only to orbital 
stability constraints but also to other conditions, for instance, 

- the spin axis (or centerline)-station vector angle, for any 
tracking station, should be between 55° and 125°. The reason 
for this constraint are -8 db and -10 db in the antenna patterns 
in the regions bounded by centerline-station vector angles of 
less than about 40° and greater than 135° . 

- ecliptic plane apogee-sun angle between 15° and 60°, decreasing 
with time. This angle is defined as the angle between the Earth 
Sun line and the projection of the geocentric vector to apogee 
onto the ecliptic plane. In other words, the projected apogee 
vector will point to the sub solar point after between 15 and 60 
days after injection. 

For the fast mapping of the launch window, these authors used 
the' Approximate method of SABAC . Diagrams such as Fig. 2.14 were 

produced* quoting from [2-9]: 

"The rapidity of this program allows one to map out a complete 
launch window in a single computer run of less than two minutes. 
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whereas use of numerical integration would require many hours. While 
this program is approximate and is not intended to be highly accurate, 
it provides an extremely useful picture of the launch window as 
a basis for more detailed study. This program was obtained for the 
IMP project office (p. 2). 

On p. 4 of the same report, commenting on the SABAC lifetime con- 
tour of Fig. 2.14: 

"Various parts of this contour were chosen as test points in 
the Encke program and it was found that the contour was fairly 
accurate until approximately March 16, 1971, near 1600 to 1800 
hours, where some complex forces apparently are beginning to com- 
bine in a manner that SABAC may not consider. As can be seen 
by the points plotted on the curve, this complex action is most 
significant around March 26 and appears to be disappearing at 
about April 10 and therefore is probably a cyclic occurrence 
related to the Sun. For this reason, if the launch is to occur 
later than about March 24, extreme care must be used. Seyeral points 

plotted on Jan. 27 show just how sensitive the lifetime is to 

ll HI 

injection time where a difference of 1 15 in injection time 

means the lifetime decreases from more than 3 years to about 
4 days! Despite the above-mentioned complexities, Fig. 1 is 
an excellent starting base for a detailed look at the launch 
window 11 . 

As a final point of interest, in Kaufman and Muhonen^ study, 
a Monte-Carlo procedure to account for the dispersion at in- 
jection should be mention here. Fig. 2.15, taken from 
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[2-9], p. 13, describes the injection covariance matrix. In view 
of the significant magnitude of the off-diagonal terms, it was 
thought to be inadequate to only examine 3-a perturbations of the 
diagonal elements of the covariance matrix. The devised Monte- 
Carlo procedure generates a set of 350 random state vectors having 
a normal distribution about the nominal, as defined by the co- 
variance matrix. These 350 random vectors are then converted into 
the SABAC input coordinates, and using the SABAC program, 350 
corresponding launch windows are generated for each launch day 
considered! (Needless to say, the cost of such a Monte-Carlo 
study would have been prohibitive if carried out by conventional 
numerical integration). The launch window lower limit differed 
by no more than + 15 minutes from that corresponding to the nominal 
state vector, and 99% of the upper limits were within 30 minutes 
of nominal. On that basis, Kaufman and Muhonen could conclude 
that raising the nominal lower limit by fifteen minutes and lowering 
the nominal upper limit by 30 minutes should avoid any problems 
caused by injection state errors. 

2.5.4 Mother-daughter mission 

A satellite mission on an orbit of large eccentricity, in which 

a subsidiary satellite ("daughter") will be separated from the main 

satellite ("mother"), is at present being planned by NASA and ESRO. In 

the mission analysis of this spacecraft, S.J. Paddack, D.P. Muhonen and 
r 2— 22 1 

G.B. Fried 1 J used the approximate criteria method and SABAC to 
generate a number of launch windows, spanning intervals of several hun- 
dred days, for various values of the argument of perigee. Fig. 2.17 is 
an example reproduced from [2-22]. 
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2.6 Auxiliary Programs 

.In this section, we shall briefly describe some auxiliary programs 
developed under this grant and to be used in support of the main 
program SABAC for the development of launch window maps. They are: 

- a program called "ECLIP" (existing in Versions 2, 3) which 
defines the relevant orbital and injection parameters, for 
given radius and velocity at injection, given inclinations 
of the velocity and satellite spin axis vectors on the 
transverse. It is often desirable that, at a nominal, "ideal" 
time, the velocity at injection (Version 2) or spin axis at 
injection (Version 3) be normal to the ecliptic. 

- a plotting program for the SABAC output, called "SABPL2" (a 
slightly modified version was written by G. Fried of NASA GSFC). 

2.6.1 Program ECLIP 

ECLIP (here, more specifically, ECLIP 2) is a program designed to 
determine, on the basis of given: radius at injection , speed at injec- 
tion , inclination on the equator , flight path angle (i.e. angle between 

V. and the transverse vector, in the direction of flight), the follow- 
mj 

ing quantities: 

1) The nominal injection time , on any given day of the year, 

which guarantees perpendicularity of the velocity vector, 

V , to the ecliptic plane (the spin axis is sometimes 
inj 

assumed to be aligned on V^^). 
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2) The range of injection time , or "launch opportunity strip", 
within which the spin axis of the satellite, at injection 
lies within +5° of the negative normal to the ecliptic 
(the injection is southwards). 

3) The solar aspect angle at the boundaries of the above"strip", 
and at the nominal injection time. 


The goemetry of the problem is shown in Fig. 2.18.. Let (X £ , Y^, Z £ ) 
be the geocentric ecliptic system, and (X , Y , Z ) the geocentric 

Ot t-l Ot 

equatorial system. The coordinate transformations are 


X 


e 

Y 


e 


Z 

e 



r 

1 

0 

\ 

, 0 

' 

— 

0 

cos e 

sin e 



0 

-sin e 

cos e 


J 



j 

• 


and 








X 

a 


i 

0 

0 


Y 

a 

= 

0 

cos e 

-sin e 


NI 

p 


0 

sin e 

COS E 

) 



The following vectors and scalars of special importance: 

a) r. ., unit vector to the point of injection (R. . = R. .*r » 

inj inj inj 

radius vector at injection) 

b) v. ., unit vector along the velocity at injection (V - V *v . . , 

inj ln J in J 

velocity at injection) 
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c) P, unit to perigee and Q unit along the positive semi-latus 
rectum 

The dependence of P on r ^ n j » ^i n j reca H ec ^ here. 

-¥ 4 - 

R. . = R. . cos Vj„. P+R . sin v , . Q 

mj inj ln J inJ xnj 


V. « = /u/p [-sin v. , P + (e + cos v. . )Q] 
inj 1 mj v mj 


Solving for Q, 
1 


•+■ -*• 
Q = r 


inj sin v ±nj 


COS V . . 

mi f 


sm v 


inj 


Then 


? inj * - -W ? 


1 + e cos v. . e + cos v 

iai + — : 3 sl ? 

sin V. . F sm v inj 


mj 


inj 


Finally, with C = | r A r * | = /pp 

sin v , 


e + cos v , . ^ 

± mi ■> 

P ^ 

1 + ecos v. . inj 


Jdil— - V. 


1 + ecos v. . p inj 
xnj 


( 2 . 6 - 1 ) 


d) y, flight path angle, is the angle between the velocity vector 
at injection and the transverse vector, in the direction of 

~y 

flight. Y is taken to be > 0 if the projection of V on R^ 
is positive , i.e. if injection occurs after perigee. 

4 - ^ 

e) Y » the satellite spin axis angle, is counted from V. . into s, 

s’ r mj 

unit vector along the spin axis (assumed to be in plane r. . , . ), 

xnj mj 

positive in the direction of positive y* 

It is assumed that the injection is southwards. Of two possible 
values of longitude on nodes in equatorial plane, the one 
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to be retained is that which delays most the appearance of 
long eclipses (more than several hours , near apogee) . As 
will be seen later, this amounts, for posigrade orbits (i £ 
and an injection in the Northern Hemisphere, to requiring that 
ft^ (which could be either in the fourth or second quadrant) be 
in the second quadrant (or Q + tt, longitude of descending node is 
in the fourth quadrant) 


Numerical injection time 

By definition, the nominal injection time is that time t ^ , at 

inj 

which the velocity at injection (or vector v^ ) , is aligned on the 

negative normal to the ecliptic. Thus, v. . = [0 0 -1) = [0 sin £ -cos e] 

ltlj E U 

Evaluating r A v .in the £-system, if n is the unit normal to the orbit, 
inj inj 


r . A v . . = cos yn 
inj inj 


Thus 


( r in j) e = [ny cos y, - cos y , -sin y] ( 


or 


( r in j ) a = ["y cos Yt “"x cos y cos e + sin y sin e, 


-n^ cos y sin z - sin y cos e] Q 


( 2 . 6 - 2 ) 


Therefore (r. . ) can be determined from y, once n is known. Given 
inj a 

i , n is related to through 
a’ a ° 

(n) = [sin ft sin i , - cos ft sin i , cos i 1 

v ' L a a a a a a 

= [sin ft sin i , - cos c cos ft sin i + sin e cos i , 
a a a a a 


cos e cos i + sin e cos ft sin i ] 
a a a e 


(2*6-3) 
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The problem is now simply to find ft such that v, .is normal to the 

a inj 

ecliptic. Let (RA) £ be the celestial longitude of the projection of r^^ 

on the (X , Y ) plane. Then 
0 0 

n = [- sin RA £ , cos RA £ , 0] £ 

= [-sin RA , cos e cos RA , sin £ cos RA_ 1 

£ £ Ot 

i being imposed, if z is the unit along the Z -axis, 
ot ot 

-y -y 

n-Z = cos l = sin e sin RA 
a a e 

The quadrant for RA £ (of two possible) being chosen from shadow considera- 
tion, as explained above, the angle is determined from 
cos i 2 

cos RA = — : — sin RA = + sqrt(l - cos^RA ) (2.6-4) 

e sin e e — e 


where i a , e are known. 

Now from Equation (2.6-3) 
n v 


cos i 


cos ft = “ 

a sin l 


a cos 


sin i sin 
a 


£ 

£ 


and ft is in the same quadrant as RA . Note that ft has been obtained from 
a £ a 

Equation (2.6-4) and (2.6-5) only, and is thus (as well as n) independent 
of y , flight path angle. The injection time, t^ , is obtained as follows. 

Let B, be the angle, measured positively eastwards, between the injection subpoint 
on the equator, and the orbital descending node (Fig. 2.18). From vectorial 
equalities, if <5 is the injection point, 

sin £ = — P — 1 cos £ = + sqrt (1 - sin 2 £) (2.6-6) 

tan l 

a 

Now let (2.6-6) be computed for the launch point (yielding ^ ) and the 
injection point (yielding £ 2 ). The Greenwich hour angle at the time of 
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injection is 

G. , - G^ nn __ _ _ „ - lljn + 15*041688* (d*24 + t, .) = G 0 + 15.041688 (in deg) 

inj 00:00 H.U.T*,on day "d" inj' 

t. . * in hours, has to be comprised between 0 and 24* (2*6-7) 

mj 

Also, if $ is the longitude of the launch site, 

L 

G, ,=$+£,- + Q + 180 (degrees) (2.6-8) 

inj L 1 ot 

t , is computed from 
inj 

‘inj ' 15?041688 (hours > £rom 00:00 hrs “• T ' ) 

15.041688 (in deg.) 

s. 

Multiples of 24 are added or subtracted to the numerator so that t^. 
lies between 0 and 24. 

The reduction to perigee, necessary for running data generated by 
ECLIP in Program SABAC, is embodied in formula (2.6-1), in which, 
successively, 

C 

e 

P 

a 


R .V. . cos y 
inj mj 


sqrttl + 2(^p 2 * (~ V inj 


_ _JL 


R. . 


)] 


U 


1 - e z 

jf 2 


= 2 


/F 
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h A = a(l + .e) - 

hp = a a - e > ” % 

2 „ 

cos v. . = (- - R. .)/(e*R, .) and sin v, , = (Sign y)sqrt (1-cos v ) 

inj y inj inj inj inj 

from (2.6-1), is equal to [X p , Y p> Z p ] in the a-system. 

Z p , longitude of perigee, is given by: sin£ p = Z p ; 

cos S p = + sqrt (1 - Z p ) 

co S E ln . - e + (E. nj cos v lnj )/a ; sin E. n . = sqrt (1 - cos 2 E. n .) 

3/2 

T . = ■%= — (e sin E, , - E , .) 

injection to perigee vp inj inj 

(This quantity is positive for injection before perigee or y < 0; 
negative for injection after perigee or y > 0) 

^ = "sin - i~ n cJ * 

^ a 

(u argument of perigee, is obtained from cos = 1^'£; sin = 

"•[V* !• 

Launch Opportunity "Strip" 

The launch opportunity strip is obtained by setting limits on the 

■V 

angular departure between the satellite spin axis vector, s, and the 
negative normal to the ecliptic. Now, assume that the nominal injection 
time has been determined as above, 6, latitude of injection, being fixed, 
at time t (t . , is supposed to be the nominal injection time), we have 
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r, , ^ = [cos <$ cos RA , cos 6 sin RA , sin 6] 
in j , t a , t a,t 

n t = [sin i a sin 0^, - sin i a cos n^, cos ij 

fi a,t = *Vt* . + 15. 041688 (t - t* ) 

’ inj mj 

RA = (RA) * + 15.041688 (t - t* ) 

a,t a ’ C inj J 

The unit vector along the spin axis of the satellite is given by 

1 = cos (rhfg) n t A r inj>t + sin (Y+Y s ) 

The (spin-axis, negative normal to the ecliptic) angle a, is given by 

sin o = |S A Z I (0 < a < 90°) (2.6-9) 

1 e 1 

and the solar aspect angle 8, given the Earth- Sun vector ^ E _g a t hy 

sin 6 = |1 E _ S A S| (2.6-10) 

An example of result of these calculations is shown in Eig. 2.18. 

A Fortran V computer program, called ECLIP, has been written at C-MU 

and is described in documentation D-2.lt determines, for given > V inj* 

i , v, v . the nominal injection time, the launch opportunity strip, the 
a s 

values of angles ct and 8. It has also been used by S.J. Paddack in the 

[2-8] 

mission analysis of IMP-G 

Modified version (ECLIP3) 

In this version, the nominal injection point is defined on the basis 
of the satellite spin axis vector s (defined by y ) being normal to the 
ecliptic. Formula (2.6-2) still holds, with y g replacing y. Obviously, 
the latitude and longitude of injection will depend on y . The time of 
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injection, computed as above, only corresponds to the time when the 

orbital plane (of fixed i^) contains the normal to the ecliptic; it 

is thus y - and y- independent. For instance, if we chose y = -y, 
s s 

s is along the transverse and the injection occurs in the ecliptic. 

With that equality, the launch window strip can be determined without 
any need to restructure program ECLIP. 

2.6.2 Program SABPL2 

SABPL2 is a plotting program documented in D-3, accepting the 
punched output from SABAC 1 or SABAC 2. It plots the launch window 
and launch opportunity "strip" as defined above. The strip is hori- 
zontal (if output is from SABAC) or oblique (if output is from SABAC 2, 
with the "strip" defined by ECLIP) 
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FIG. 2-1. GEOMETRY AND NOTATIONS. 





HOUR OF LAUNCH 


INITIAL ORBIT: 259.000 km; h p : 250 km; 



1 11 21 1 11 21 31 10 20 

DATE OF LAUNCH 


FIG. 2-3. EXAMPLE OF LAUNCH WINDOW MAP 
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(x 10 8 ) 



TIME FROM LAUNCH (DAYS) 


FIG. 2a 



FIG. 2< 



FIG. 2b FIG * 2d 



FIG. 2e TIME FROM LAUNCH (DAYS) 


FIG. 2f 


FIG. 2-4. VARIATION OF ORBITAL PARAMETERS WITH TIME 
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0 100 200 300 400 

TIME FROM LAUNCH (DAYS) 


FIG, 2'r5» Short-term and long-term effects on altitude of 
perigee. 



SEMI-MAJOR 



FIG. 2-6. CONSTANCY OF SEMI-MAJOR AXIS. 
(ORBIT OF FIG. 2-4) 



TIME FROM LAUNCH 


* 

FIG. 2-7. SHORT-TERM AND LONG-TERM CHANGES OF e AND h. 






FIG. 2-8. LIMIT ON AMPLITUDE OF MODULATION 



ECCENTRICITY 



FIG. 2-9. VERY-LONG-TERM EVOLUTION OF THE ECCENTRICITY. 
(After B. Shute' ^ . 












FIG. 2-12. COMPARISON OF RESULTS FROM NUMERICAL INTEGRATION AND 
THIS THEORY. ' 
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Nominal IMP-1 Launch Window (0 km Initial Perigee Drop) 

FIG, 2-15 (From Ref. [2-9]. 



INJECTION STATE VECTOR 


1 . ALTITUDE, ABOVE MEAN EQUATOR (km) 

2. VELOCITY (km/sec) 

3. FLIGHT PATH ANGLE (deg.) 

4. FLIGHT PATH AZIMUTH (deg.) 

5. LATITUDE (deg.) 

6. LONGITUDE (deg.) 


231 .35075 
10.806989 

0 . 

73.5944 

-23.966 (SOUTH) 

1 1 1 .449 (EAST) 


COVARIANCE MATRIX 

© © •© 

319.475 -.301558 -1.89925 

.000314949 .00168757 

.0536598 


@ 

© 

© 

3 .06964 

-1 .85248 

-7.01463 

-.00303874 

.00174996 

.00661598 

-.0191042 

.0116318 

.0393667 

.0982999 

-.0183223 

-.0649491 


.0111156 

.0392638 


.159351 


FIR. 2-16. COVARIANCE TABLE, IMP-I ORBIT 
(FROM REF. [2-9] ) 



KOTHER/DAUGHTER RUSSIAN OMEGA = 0 7/77-2/73 



DAYS OF THE 1972 (REF. JAN. 1=0) 


FIG. 2-17. ■ LAUNCH WINDOW OF MOTHER-DAUGHTER MISSION (From Ref. [2.22]), 
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CHAPTER 3 

A Study of Orbits of Large Eccentricity Quasi-Normal 

to the Ecliptic 


3.1 Introduction 


Relatively early in the present study (1968-1970), it appeared of 

interest to initiate a study of the "practical" stability (in the 

sense specified in Chapter 2) of high eccentricity orbits having an 

inclination on the ecliptic close to 90°. This was to be the case 

for satellite IMP-G, which was launched in June 1969. A very detailed 

[ 3 — 1 ] [ 3 - 2 ] 

description of the results has appeared elsewhere ’ 

A thorough study of IMP-G orbit and launch time was carried out at 


r 3— 3 1 

NASA GSFC by S.J. Paddack 1 1 , who used ascomputing technique ITEM 
and a Perturbation Routine for final, accurate results by numerical 
integration on a computer, and our program SABAC, based on the 
approximate criteria approach, for the fast generation of global launch 
windows. ECLIP, written at C-MU, was also used to define the launch 
opportunity strip, defining a nominal time, and an interval on both sides 
of this nominal time, in which the alignment of the satellite spin axis 
on the normal to the ecliptic is closely realized. 

As is schematized on Fig. 3.1, taken from Ref. [3-3], it is clear 
that if, to simplify the reasoning, we assume that the velocity at injec- 
tion is very nearly coincident with the satellite spin axis and is normal 
to the ecliptic, the resulting orbit will be inclined by i = 90° on the 
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plane of the ecliptic. The corresponding inclination on the equator 

i , for a southwards injection, will be comprised between 90° (if 

Q = Q = +90°) and 66.55° = 90° - e (if U = U = +180°). The fourth 
e a e a 

and first quadrants for ,£2 would correspond to retrograde orbits, 
and the third one might lead to prolonged times in eclipse, near apogee, 

r 3 — 3 1 

for up to 9 hours 1 J : therefore, these three quadrants are not to be 
considered, and we shall assume from now on that 90° i or < 180°. 

In the following, the approximate criteria shall be used consistent- 
ly; - the lifetime criterion is described in Section 2-4. 

[3-1] 

3.2 Study of the Stability of Orbits Nearly Normal to the Ecliptic 


3.2.1 Simplified model: planar case. 

As an approximation, we shall first consider, that in view of the 
smallness of the inclination of the moon's orbital plane on the ecliptic 
(I^= 5.145°), those two planes are approximately coincident. With the 
notations of Chapter 2, let i be the inclination of the satellite orbit on 
the orbital plane of perturbing body "d". As said above, we have 
approximately (Fig. 3.2) 

in = i g - 90° 

Long-Range and Very Long-Range Stability 


As explained in Chapter 2, we shall use Lidov's "11", secular theory. 
With the above simplification, we can define "A^" as the sum A^ + Ag 
(Note again that A^/Ag is independent of a, and equal to 2.18). 

The changes in the orbital elements due to the sun and the moon, per 
satellite orbit, are (angles, referred to the ecliptic, are not subscripted) 



3-3 


6a = 6i = 6ft = 0 

6e =ee 1/2 (A M + A g ) 
^ A S . 


sin 2to 


6to = 


20 


4 = (5 cos 2to - 1) (1 - 


') 


5 cos 2to - 1 


c 2 = (1 ~ e o )(5 cos 2co o - 1) 


( 3 . 2 - 1 ) 

(3.2-2) 

(3.2-3) 

(3.2-4) 


Subscript n 0 M refers to initial values* 

To this approximation, orbits initially normal to the ecliptic will 
remain so for all time and will have a constant longitude of nodes* Now 
consider Lidov's constants, c^ and C 2 , which are integrals of the 11 , 
secular differential equations of motion: 

Cl = cos 2 i = e a cos 2 ± o = 0 |f 

_ ■ - c 2 
c 2 = (1 - e)(| - sin 2 to) = (1 - 0(f - Sin 2 u 0 ) = (3.2-5) 


Any of the orbits will be represented by point (c 2 , 0) on segment AC of 
the c 2 axis of Lidov's (c^Cg) diagram (Fig. 3.3). From Equations 
(3.2-1) to (3.2-4), the evolutions of to and e are described by Lidov's 
discussion ^ 3-4 ^ , in which A d is replaced by A M + A g . 

Let <o* = i arc cos(|) = 39.23°, = tt - io 1 , co 3 = it + Uj, to 4 = 2 tt - 


Then, considering that 6e = -A(l - e)e 


I/? sin 2oy 


A 
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a) to^ < ui < oj 1 (Region 1, Fig* 3*4) or 0)2 < w < CO3 (Region 3, Fig* 3-4) 

If initially, w 0 < 0 (os 0 <71), e decreases until to reaches 0 (tt). 

The minimum for e is given by 


1 2 
- — 


e . - 7 er (5 cos 2to - 1) 

min 4 0 o 


(3*2-6) 


Thereafter, e increases and reaches e at to . With R , 

7 max e,max # 

& 

earth's radius, and h^, critical height of perigee, it is obtained 


from 


R~ + h 


1 - e = 1 - [1 - 

max 


(3.2-7) 


5 cos 2co - 1 = (1 - e )(5 cos 2co - 1)— 5 (3.2-8) 

e ' max * <4* 


and co is in the same region as to . If initially, (0 o > 0 

e,max 0 

(to o > it), there is never a decrease in eccentricity. From the 
viewpoint of long-range stability ( Criterion 1 of approximate 
criteria method), sub-regions la, 3^ are acceptable. 

b) <0* > co > co* (Region 2, Fig. 3.4) or co* < co < co 4 (Region 4, Fig. 3.4) 

ir 3*n" 

Here co decreases. If initially, co 0 > ^ (co > Tp- ) , e decreases 

TT 

until to = y (3tt/ 2). The minimum for e is 

e . = ^e 2 (l - 5 cos 2to ) (3*2-9) 

Thereafter, e increases up to the value given by Equation (3*2-7). 

If initially w 0 < (co 0 < 3tt/ 2), the eccentricity is always in- 
creasing. Thus, from the viewpoint of long-range-stability, region 
2b. and 4b. are acceptable. 


/ 
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By reason of symmetry, and since, as was mentioned earlier, our 
interest here lies in orbits with southwards injection, the analysis will 
be restricted, without loss of generality, to the second and third quad- 
rants of the orbital plane (Fig. 3.2). On Fig. 3.4, point "B" corresponds 

to to = (02 and appears as an unstable point, whereas point "D", for which 
& 

to = to 3 , appears as a stable point. Only 2b and 3a are acceptable for long 
range stability in these two quadrants. 

If very- long-range stability is now considered, an assessment of the 
orbital lifetime can be made here, provided it is fairly large compared to 
the periods of the perturbing bodies (in practice, the required lifetime 
equals many orbital periods of the moon, but only one or a few orbital 
periods of the sun.) This so-called "long-range" (LR) lifetime reads, 
in satellite periods. 


Region 2b (ir/2 < m 0 < m 2 ) 


(-) 

'Vlr 


20 ,*/ 2 . 
T~ f C 1 dm 


m„ 


( 3 . 2 - 10 ) 


Region 3a (m 2 < m 0 < r) 


(-> 

Vlr 


20 


P c _1 dm 


( 3 . 2 - 11 ) 


in which £ is given by Equation (3.2-4). As an example, L is given, 

in days, as a function of w, for an orbit of e c = 0.945991, a = 124,283 km, 

x = 5.0468 days (Region 3a., neighborhood of the ecliptic) (Fig. 3.5). 
sat 

* 

If the orbit originates at co= m 2 + n (n small and positive angle) 

evolution B -*■ A will lead to a larger L TT , than evolution B,-> C: 

— LK T 
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if X is the point having argument ir - 2to i = 101.54°, the time spent along 
B_ X is the same as that along B + C. Thus 

t B-^A = t B->-X + fc X->A > t B^C 

In conclusion, larger lifetimes will be possible in region 2b, the larger 

* _ _ 

the closer the argument of perigee is to In region 3a, ror maximum 

a) should be made equal to o >2 + 0. The upper limit for in this 

* 

region is infinite since L^-*- » when u •> W 2 * 

SHORT-RANGE AND INTERMEDIATE-RANGE STABILITY 

The short-range behavior of the eccentricity, which determines the 
short-range stability, is described by Equation (2.4“7 ) 


6e SR = 


x /2 r 


A, 


• [ “^77 3 3,M + ~V^ 6 3,S ] 

E M 6 S 


(3.2-12) 


From Fig. 3-1, it is apparent that if 10 were 180°, as is the case for an 
injection, at perigee, in the ecliptic plane, E ^ wou ^-^ vanish and ini- 
tially, the short-term stability would be neutral (Criterion 2). 


Due to the short-range increase of to, as given by 

\ , P M 3 A q P<? 3 

— Tfz + { 

e m m e M ' " *s 


§> 3+ <r4fi>§) 3 ] >° W.2-13) 


there is, however, intermediate-term instability (Criterion 4) since £ 2 ^ 
will become < 0 at the next orbit. Thus, in the quadrants considered, 
for short-term stability it is required that the perigee be above the 
ecliptic plane 
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2 < " < * 

which also ensures LR stability (Criterion 1) 
1 1/2 

ie Tn = 7 ee A sin 2m < 0 
LR 4 


(3.2-14) 


(3.2-15) 


For intermediate-term stability, a margin (of the order of a few degrees) 

should be provided, so that the perigee is sufficiently above the plane 

of the disturbing bodies. In order to obtain actual lifetimes which are 

as long as possible, one could possibly select those launch days in the 

year leading to a slope, on the e vs. time from launch curve, and over a 

time of the order of x /4, which is assmall as possible. For fixed e D , 

s 

ij = m 0 and thus fixed values of Lidov's constants c* and c 2 , (^e^ R ) M 
is fixed. One requires to make < 6 e gR g > Tg ^ as small as possible. As 
an example, for a celestial longitude of the radius vector at injection 
RA £ , in the fourth quadrant (Fig. 3-6), the best launch day of the year 
will be that for which the unit vector to the sun is along axis 3 r = 3^ 
normal to the orbit, on the average over x s /4 (half a solar modulation.) 
An example is treated in Section 3.3. 


Conclusions from simplified model 

In summary, to the approximation of the simplified model, and with 
an analysis restricted, without loss of generality, to the second and 
third quadrants of the orbital plane, it is concluded that 

a) the highest realizable value, in region 3a. of Fig. 3.4, for 

the long-range assessed lifetime is defined by Equation (3.2-11), 



3-8 


where io 0 is the minimum feasible to. 

b) short-term and long-term stability are strictly realized for to 
in the second quadrant, whereas intermediate-term stability 
requires that to differ from ir(or it/2) by a negative (positive) 
margin, in practice a few degrees for a lifetime of one year. 

c) it is possible to define a best day in the year leading to maxi- 
mum actual lifetime for given e 0 ,w and inclination on the equator. 


3.2.2 Effect of the inclination of the moon's orbit on the ecliptic 


We still assume that the satellite orbital plane is normal to the 

ecliptic, RA is defined as the celestial longitude of perigee (Fig. 2-11) 
e jP 

If P is sensibly in the plane of the ecliptic, and above the moon’s plane, 

let £2 be the longitude of nodes of the moon's orbital plane, referred to 
e,M 

the ecliptic. Thus 

e,M e,P e,M 


Stability is assured in the long-range, since ( 5e ) LR> s " 0 and 

< o. The locus of southernmost (northernmost) admissible perigee 

v 'LR,M 

points on the unit sphere is obtained by writing 

(6e) LR fc (1 - e)e 1/2 (A^ + Ag* g ) = 0 (3.2-16) 

with Til = TT - u,. It is approximately the arc of great circle having 
d d 


normal Z„ inclined by 

A m 


i = 




+ K 


hi 


(3.2-17) 


on the normal to the ecliptic, Z^, in plane (Z^,, Z^).In the long-range, 
therefore, 1^ may be accounted for by rotating the ecliptic by i about the 
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nodal line of the moon. This new plane of reference, called (o~, then re- 
places the ecliptic. 

As an example, for a nominal injection at perigee and in the ecliptic 
(fourth quadrant), Fig, 2.11 shows that an inclination on the equator 
i Q = 90° should maximize the long-range assessed lifetime if = 0° 

(which is in the case in early spring 1969), in a range for i^ 


Z. 

2 


c < 




An example is treated in Section 3.3. 

In the short-range, were the sun alone, condition 
IT ^ 

Y < W < IT 

would still hold, i.e. stability in the short-range would be realized 
when the perigee is in the second quadrant of the orbit. If the injection 
occurs at perigee, in the ecliptic plane, the short-term effect due to 
the sun alone is zero; the moon critically determines short-range stabil- 
ity. If the projection on the ecliptic of the vector to perigee, OP, in 
Fig. 2.11, is normal to the moon’s nodal line in the ecliptic, ON^, the 
moon is certainly favorable or neutral if 


£ E > 
*1,M 2,M 


0 


and cannot be satisfied throughout the lunar month. 

In the intermediate-range, the margin on m, to which we referred 

above, should not be construed with reference to plane w.. . The condition 

2 

on the best day of launch still holds approximately, to 0(1^). 
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3.2.3 Effect of a small departure of the orbital inclination from normal 
to the ecliptic 

Earlier, we defined a "nominal" orbit as one normal to the ecliptic, 
and it is desirable to qualify the effects of a slight departure, Ai £ , from 
being normal to the ecliptic. Such a departure will be caused, for 
example, by the Earth’s rotation for a launch slightly earlier or later 
than nominal. 

r 3—4 1 

Lidov's formulae, in the "11", secular theory 1 , written up to 
0(Ai £ ), yield 

6e LR = ~\ A(1 " e)el/2 - S ~ V ~ (3.2-17) 

For long-range stability, it is required that Se > 0, or sin 2w < 0. 

LK 

To the. same approximation , 

Srn = ^ Ae ^ 2 [i - sin 2 a>] (3.2-18) 

LK L j 

Therefore, the developments of the two previous sections concerning long- 
range stability apply. 

In the long-range, the orbital inclination varies according to 
1 Ai 

Si =-jA(l - e) sin 2 w — 17 — (3.2-19) 

LK 4 1 / 2 

7T TT 

and i will increase (decrease) for Ai £ > 0 or i < ^ (Ai £ < 0 or i > 7 ) 

7 r 

and tend to -j for stable orbits. 

The rotation of the line of nodes will be of order Ai £ , 

6Q LR = - \ A Ai e t(l ~ e)sin 2 w + |] --77 (3.2-20) 

Developments relating to short- and intermediate-term stability involve 

2 

geometrical conditions in the orbital plane, and still apply to 0(Ai £ ). 
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з. 2.4 A Priori Prediction of Orbital Lifetime 

It has been seen in Section 2. 4. 3. 3 that, in order to assess the 
orbital lifetime, T^, and compare it to the required lifetime, L, one 
possible method consists in computing Lidov's constants c^^ ^ and c^ 
for the sun and the moon, and to weigh them with amplitudes A^ and A g 
to obtain resultant c^ and C 2 • It has also been mentioned that for 
orbits nearly normal to the ecliptic, i^ ^ ig ^ 90°, the approximation 

и, , ^ co„ might not hold at all if, say, tt - cu, is a small angle (injec- 

M S a 

tion southwards, near the ecliptic). In particular, in the equation 

1 . 

6e r „ J = A j eG /Zsin2i J SiT1 2aJ J 

LR,d d d d 

sin 2oj d might be of different sign for the sun and for the moon. 

Therefore, it appeared necessary, for satellites spending a signi- 
ficant fraction of their lifetime "between" the orbital planes of these 
two perturbing bodies, to come up with a better method of estimating the 
orbital lifetime T A . Typically, when the unfavorable influence of the 
sun (6e > 0) on the evolution of the eccentricity, as would be the case 

if c i , C 2 were with reference to the moon's orbital plane only, it was 
(3 2-1) 

found that the inaccuracy in some parts of the contour of the map 

(determined by the lifetime criterion) was of the order of 0.3h, a 
large error in the case of IMP-G. 

The alternative method adopted in this case has been briefly de- 
scribed in Section 2. 4. 3. 3, but will be repeated here. Plane as defined 
above, is the reference plane used to compute the argument of perigee of 
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of the satellite uu (Fig. 2.11). Now, half of the orbital lifetime, 
T^/2, will correspond to the time needed for the perigee point to reach 
plane . The time-rate of change of is computed in equation 


_ <Sco-w dm 

Jt~ = SsT” ^LR,M + ( 6 t 'LR,S 


(1 + w) 


as being due to the long-range effect of the moon and, in order to cap- 
ture the linear trend of vs. time, the average of the short-range effect 
of the sun Or^)™ 0 taken over half a solar modulation cycle, t s /4 (see 
Fig. 2.5), As an example, it is interesting to note, in Fig. 3.7, where 

o) vs. time has been obtained from numerical integration program EOLA, 
a 

that the crossing of plane 5k (as marked by the arrow) quite accurately 
corresponds to the topping off of the height of perigee. 

To be more specific, we consider the "11”, short-term theory, as 
embodied in Equations (2.3-14). Per orbit of the satellite, 5T sat > if e s 
is neglected (a^ p g ) and i g % | , for body "S", 


6co„ 


dr 


sat 


1/2 

e «„s 


In the vicinity of the ecliptic, 


2 



2 

£ _ can be neglected. 

2 , S 


Now 


2 

5 

1,S 


Si cos 2 (RA g - RAp) 


in which RA p is the celestial longitude of the perigee. Let <j> = RA g - RAp. 
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,6to- 


v r v» _ A 

6t r g /4 — 5 S 
sat s 


* = A„ e ' <4 


1/2 , 1 + cos 2<f> 


to/4 


- 5 A s e ' /2 + \ A S £ ’ /2 x <- I sin 


- > + ^ (A) ~ \ w 

1 6t ,/LR,S V 6t /LR,S 
sat * sat 


in which 

w 


4 . 

- sxn 2<p 
ir 


inj 


For an injection in the neighborhood of the ecliptic, and given i a , a "best" 

6 GO 

day is one which minimises or for which 


(RA g - RAp)^ - $ = ± k ~ (k positive integer or zero). 


An example is given in the Section 3.3. In the same section, Table 3.3 I 
also shows the good agreement obtained, by this procedure, between the 
predicted lifetimes, the latter being obtained from numerical integration 
programs (EOLA or NASA's ITEM). 


3.2.5 Effect of the Earth's oblateness on the orbital lifetime 

In the course of the present study and the subsequence application to 
IMP-G, it was found that orbital lifetimes can be significantly enhanced 
by the effect of the equatorial bulge (J 2Q term in the Earth's potential), 
by up to 20% in some cases under study. It is recalled that so far the 
Earth's potential had been considered spherical in the analysis. As is 
well known, due to J 2Q , there (till be no secular changes of the satellite 
semi-major axis, inclination or eccentricity. The line of apsides (in the 
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orbital plane) and the line of nodes (in the plane of the equator) will 
rotate at rates proportional to 4-5 sin 2 !^ and cos i^, respectively. 

Since all i^ considered here are higher than critical, is < 0, 

and in the range Investigated, its magnitude is maximum when i^ = 90°. 

The same "beneficial" effect of oblateness on the stability of high 
eccentricity orbits of natural satellites is mentioned in examples given 
by J. Kovalevsky ^- 3 ^ and Lidov ^ 3 5 ^. 

Typically, the kind of orbits studied here have perigees which rise 
very little (10 3 to 2 x 10 3 km) over the whole orbital lifetime. This 

is in contrast with the more frequent occurrence described for example 

, [ 3-6 1 

by Shute 1 4 . 

3.2.6 Conclusions of the study 

The conclusions and practical implications of the above study, when 
applied (without loss of generality) to a satellite launched southwards , 
into an orbit quasi-normal to the ecliptic, with a perigee in region 3a 
of Fig. 3.4, are as follows: 

a) High celestial latitudes of the perigee are required for the 
stability in all ranges. They will be the more favorable the 
closer the argument of perigee referred to u)„ is to ir - j arc cos( 
In particular, a positive flight path angle (i.e. injection after 
perigee) will be beneficial, within limits allowed on the drop in 
perigee height as compared to injection height, and mandatory if 
the injection is to take place in the close vicinity of the ecliptic. 
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b) For a nominal launch, at perigee and in the ecliptic, and fixed 
eccentricity e and inclination i a , it is possible to define a 
best day in the year giving the longest lifetime. 

c) The most suitable inclination on the equator, for a nominal 
launch, is that corresponding to an orbital ascending node at 

tt 

+ ^ from the moon's node. 

d) If the nodal line of the moon is sensibly aligned on the vernal 

line (the case in early spring 1969), the angular height of 

perigee above is tt - w + ip (Fig. 2.11). Therefore, if 

e,r 

< 90°, advantage can be taken of the Earth’s rotation to 
increase this angle, and consequently the lifetime, by launch- 
ing earlier than the nominal time. 

These conclusions were used with profit in the mission analysis of 

a high, eccentricity satellite in an orbit nearly normal to the ecliptic, 

r 3_2 1 

IMP-G . This study is described hereunder. 

3.3 Application to an Actual Satellite: IMP-G 
3.3.1 IMP-G orbital data 

The abovementioned study of IMP-G orbit and launch time, carried out 

r 3 3 1 

by S.J. Paddack , should be referred to for more specific details and 

mission analysis studies. Our motivation here was to use the results of 
Section 3.2 and apply them to satellite IMP-G, in order to possibly pre- 
dict the qualitative and quantitative effects of the launch parameters 
on the orbital evolution, and more specifically the orbital lifetime (re- 
quired to be larger than 1 year, even for the 3a velocity dispersion orbit). 
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A description of a typical orbit follows: 

IMP-G 

h. = 235,463 km 

h = 343 km 

P 

e = 0.946 

T sat = 5 ’ 05 dayS 

Y (flight path angle = angle (ly positive if 

% > 0): -2° to +2°. 

v in j inj 

Y g (satellite centerline angle: angle ^ spin ax i s » t0 +2 °* 

i^ (inclination of orbit on ecliptic): about 90°, at nominal time 
Injection in 4th quadrant of ecliptic 

The above list calls for a few comments: IMP-G is spin-stabilized, 

without active attitude control. It was desirable that the spin axis 

vector t , aligned within a few degrees on the velocity vector at injection 
s 

V. . , be normal to the plane of the ecliptic, within a narrow tolerance 
inj 

A£ = + 5°. Injection is made very close to perigee (within a few degrees). 

Hence, the resulting orbit will be very nearly normal to the ecliptic. For 

example, if y - y = Ac = 0°, i.e. for an injection at perigee with velocity 

and spin axis vectors exactly aligned on the negative normal to the ecliptic, 

-Z^ , the perigee at the so called nominal time will be in the ecliptic, at 

celestial longitude RA . . depending on the inclination on the equator, 

e , inj 

i (Fig. 3.8). For obvious reasons, posigrade orbits are preferred (i < 90°), 
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and to avoid long periods in the earth's shadow at the outset of the 
mission, only injection in the fourth quadrant of the ecliptic are to 
be considered: 


-90° 


< RA 


e,inj 


< 0 ° 


The constraint on the alignment of the satellite centerline on -Z limits 
the launch opportunity to a "strip" of width equal to about 1.7 hour, 
symmetric about the line of nominal launches (Fig. 3.9). As illustrated 
in Fig. 3.10, the degrees of freedom in choosing a "suitable" orbit, with 
a special emphasis put on achieving larger lifetimes, are 

a) hour of launch, HL, inside the strip, on a given day 

b) day of launch, DL 

c) inclination on the equator, i 

d) flight path angle at injection, y 

e) satellite centerline-velocity vector angle, y g 

The conclusions of Section 3. 2. 3. 6 will now be used in a systematic in- 
vestigation of the effect of these parameters on the orbital evolution. 


3.3.2 Parametric study of IMP-G 
3.3. 2.1 Launch opportunity strip 

Above described ECLIP program was used to define the launch opportunity 
strip based on a specified maximum angle A? between the spin axis and the 
normal to the ecliptic. This strip defines a range of permissible injec- 
tion within the specified tolerance. The "backbone" of this strip is the 
time of nominal injection times, at which the spin axis and the normal to 
the ecliptic are exactly aligned (Fig. 3.9). This can correspond to a 
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nominal injection, i.e. in the ecliptic plane. The misalignment Y g 
between the velocity vector at injection and the spin axis, should then 
be compensated for by an equal and opposite flight path angle, y (Fig. 3.10). 

3. 3. 2. 2 Launch window 

The stability analysis of the orbits derived from program ECLIP was 
carried out by means of program SABAC. 

a) Influence of the flight path at injection, y 

To an increase in the magnitude of y for constant R ^ n j » ^i n j * 
there corresponds a drop in the height of perigee equal to about 
4 km/deg, of change in y, in the range |yj — 2 . The rate in- 

creases with increasing y. 

All things being equal, a positive flight path angle (injec- 
tion "after" perigee) causes a high initial angle of perigee 
above plane uL , consequently a larger lifetime. This leads to 
an improvement in the "quality" of the launch window, as 
measured by the area covered by the "success" region within the 
launch opportunity strip. Fig. 3-11 to 3.12 graphically portray 
this for the IMP-G satellite. 

b) Influence of launch time on a given day 

It is apparent for Fig. 3.10 to 3.12 that the launch window 
seems to be more favorable at times earlier than that of the 
nominal injection. Fig. 3-13, which is a plot of the predicted 
lifetime as a function of the time of day, also indicates the 
same effect. As was mentioned in Section 3.2, in early spring 1969, 
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the nodal line of the Moon is very sensibly aligned on the 
positive vernal line, axis X £ (in Fig. 2.11). The angular 
height of perigee, approximately equal to it - w + , is 

t 9 r 

increased, for given i , when the injection occurs earlier 
than at the nominal time, due to the rotation of the Earth 
in intertial space. It is also clear that if i = 90° and 
for this position of the Moon’s nodal line, the lifetime will 
top off at nominal injection time, a fact illustrated by 
Fig. 3.13. 

c) Improvement of the lifetime with the satellite centerline 
misalignment angle, y s 

Figure 3.14 illustrates that the lifetime increases with uu , 
argument of perigee relative to plane . The figure is a 
plot for a nominal injection of IMP-G, and i = 90°, the high 
values of u>„being attained by the use of a negative y and a 
compensating, positive y. 

d) Influence of the launch day 

For an injection in the close vicinity of the ecliptic, a 
"best day" for given i is one for which (RA - RA ) = <j>, at 

Ot £ £ jt 

injection, is , since is minimized. This is shown by 

f 3-3 1 

Fig. 3.15, for i - 83.8° (nominal injection of IMP-G L J ). 

e) Influence of the inclination on the equator, i^ 

In the period spring-summer 1969, the celestial longitude 
of nodes of the Moon is close to 0° (within 5° over March-July 1969) . 
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Hence u)_ , and the lifetime, increase with inclination in 

the interval 66.55° £ i £ 90°. For nominal injection con- 

ditions, and w = 0, the lifetime is plotted vs. inclination 

i in Fig. 3.16. 
a 

f ) + 3o orbits 

Due to the dispersion on the actual values of the velocity 
at injection, it is important that the launch window also be 
determined for extreme cases, such as a 3c? error on the velo- 
city at injection. The probability of having more than a 3a 
error on the injection conditions is only 0.26%. The following 
list summarizes the la (1 standard deviation) with the Delta 
launch vehicle, as taken from Ref. [3-3]. 

1-a Vehicle Errors 


Latitude 

Longitude 

Altitude 

Speed 

Azimuth flight path angle 
Elevation flight path angle 
Spin axis azimuth angle 
Spin axis elevation angle 


+ 0.4337° 

+ 0.2335° 

+ 15.426 km 
+ 0.010998 km/sec 
+ 0.6526° 

+ 0.5208° 

+ 2.0435° 

+ 1.6827° 


The 3a dispersion limits on the velocity at injection were 
studied for the IMP-G launch window, and are illustrated in Fig. 3.17. 
It is seen that until day 160, approximately, the 3a window is 
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totally closed, and that around the middle of the year (from 
day 180 to 210), the window is favorable, even in the +3o 
dispersion case. 


3. 3. 2. 3 Comparison with numerical integration programs 

In more detail than Fig. 3.7, Figures 3.18 and 3.19 show the evolution 
of some orbital parameters for a sample IMP-G orbit , as obtained from 
digital integration program EOLA. The simultaneous topping off of and 
of the height of perigee have already been mentioned, and this provides an 
experimental justification to the procedure adopted to assess the orbital 
lifetime. 

Comparative lifetime values for IMP-G, as obtained from SABAC, Version B, 
on one hand, and from NASA's ITEM (Encke’s method) and EOLA (Variation of 
parameters) are tabulated in Table 3.3-1. 


DAY 

1969 

Inj . hour 
U.T. 

Y 

deg. 

V 

de g- 

Life days 
pred. act. 

Int. 

a 

prog. 

06/01 

9.488 

1.5 

-1.3 

413 

410 

VP 

06/01 

9.988 

1.5 

-1.3 

340 

370 

VP 

06/14 

9.321 

-1.28 

-1.3 

425 

404 

VP 






397 

EM 

05/01 

11.363 

1.5 

0 

362 

389 

VP 

05/01 

12.263 

1.5 

0 

319 b 

348 

VP 

05/08 

10.097 

1.5 

0 

FAIL 

364 

VP 

05/08 

10.297 

1.5 

0 

355 

369 

VP 


^P: Method of variation of parameters; EM: Encke's method 

b At T + 0.1h:350 


Table 3.3-1 Comparison of predicted and 
actual lifetimes 
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As the table shows, the average error between the predicted and 
actual lifetime values is of the order of 5%, and on the pessimistic 
side. Case 4 is a case in point. It corresponds to the evolution por- 
trayed in Figures 3.18 and 3.19. The predicted lifetime was 362 days 
and actual one 389 days. Case 7, on the other hand, illustrates an in- 
accuracy in the definition of the lifetime boundary. However, the life- 
time predicted for the next point on the same launch day (for a step of 
0.1 hour) is in good agreement with t-he actual value. 

3.3.3 Implications for IMP-G orbit 

On the basis of the above study, recommendations could be made re- 
garding the choice of an orbit having a long lifetime, ample launch 
opportunities and still fulfilling a set of additional constraints. The 
finally selected orbit would have to consider, of course, the capabilities 
and limitations of the launch vehicle. Of particular interest here, is 
the combination of a positive flight path angle with a negative "spin axis 
centerline-velocity vector" angle. The latter combination will enhance 
stability throughout the launch opportunity "strip" and/or permit injec- 
tion at more moderate southern geographic latitudes. 

3.3.4 Conclusions 

In this example, it has been shown that the method of approximate 
stability criteria could be used with profit in a parametric study of the 
influence of various orbital elements on orbital lifetimes and the launch 
window map. The analysis resulted in practical recommendations which can 
be assessed within the perspective of the global mission. 
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FIG. 3-7. ORBITAL ELEMENTS VS. TIME, FROM NUMERICAL INTEGRATION. 
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CHAPTER 4 

A Modified Lidov *s Method by Non-Numeric 
Computation, with Applications 

In this chapter we shall deal with the application of non-numeric 
computation to the theory of geocentric orbits of large eccentricity. 

The developments constitute a modification and an extension of Lidov's 
theory ^ 4-1 ^ : a modification because it recognizes that Lidov only intro- 
duces and studies five orbital elements, whereas rigorously six of them 
should appear; an extension, because the developments are pushed to 
higher order than the "11", "21", and "12" terms of Lidov, thanks to 
the labor-saving and error-free features of non-numeric computation. 
Finally, the effect of oblateness is considered and numerical examples 
are given to illustrate the degree of accuracy and the marked economy 
in computer time obtained by using the present approach. The main lines 
of the developments add a few specific examples are given here. For 
much more detailed information, the reader should refer to the Ph.D. 
thesis ^ ^ written by one of the authors of this report (R. Sridharan) , 
as the Principal Investigator’s (Marc L. Renard) advisee. 

4.1 Motivation 

It is well realized by any mission analyst that repetitive, high- 
accuracy calculations amounting in one way or another to direct, numeric— 
cal integration of the orbit, can be a very expensive proposition in terms 
of effort and computer time. Figures of the order of 10 to 15 min per 



launch "point" (day. in the year, hour in the day) of the launch window 


are quite realistic if a modified Encke's method is used on a large IBM 
system (7094 or 360 series). The awareness of this problem led to the 
development and use of a method of lower, but sufficient, accuracy, and 
of very high computational economy, which has been described in Chapter 2. 

Between these ends of the spectrum, at the inception of the present 
study, there seemed to be a real need for a theory of intermediate com- 
plexity which basically 

- would be less costly than digital integration, by a factor 
of 10 to 100. 

- would include a sixth element, thereby resulting in an improve- 
ment of the prediction of the "time flow" along the orbit 

- would be adequate for eccentricities up to 0.95 

- could be implemented on a digital computer, for the repetitive 
calculations called for in mission analysis. 

4.2 Main Features of the Approach 

There are two main features in the approach: an extension of Lidov's 

theory so as to include a sixth osculating element, and the use of sym- 
bolic manipulation on the computer. These two points are now discussed 
in more detail. 

4.2.1 Extension of Lidov’s theory 

Lidov's theory ^ 4 ^ has been described at length in Chapter 2. It is 
recalled that the equations of motion which are retained are five in 
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number. They describe the time-rate of changes of elements, a, e, i, <d, S2 
as equal to expressions in the right-hand side, in which the small per- 
turbing forces due to a third body appear as factors. But any sixth 
element, such as T A , the time of perigee at the epoch (perigee), in 
Kepler's equation of time 

E - e sin E = n(t - T^) 

is conspicuously absent from Lidov's theory. Furthermore, the ”11", "21", 
"31" terms given by Lidov for f ive elements do not suffice to obtain high 
accuracy in the case of high eccentricity orbits (0.9 ^ e < 0.95). 

Roth 33 ^developed more such terms by hand computation, but did not go 
far enough in the Legendre Polynomial (LP) expansion of the force, as will 
be shown later. Furthermore, his choice of the sixth element is apparently 
inconsistent, as will be shown in a later discussion. 

The goal will thus be to obtain a "theory" describing the perigee-to- 
perigee variations of the orbital elements. The developments will be 
rendered more accurate both by the inclusion of a sixth osculating element, 
which will result in a more accurate timing of the occurrence of perigees 
and in a better computation of lifetimes; and by carrying the LP and Taylor 
series of the perturbing forces to the order deemed necessary from esti- 
mates of accuracy. 

Alternate, relevant approaches and methods are examined by Sridharan in 
a literature survey^ - and include contributions by Kozai^ 

Musen^ 4- ** ^ ^ f Smith^ 4 ^ Fisher and Murphy^" 4 , Fisher^ and 

ri»-io] 

Cook and Scott 1 
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4.2.2 Symbolic manipulation (non-numeric computation) on the computer 

A rough assessment of the algebra and "bookkeeping" involved in pro- 
ceeding to develop a modified Lidov's theory along the lines described 
in Chapter 2, rapidly points up the need for mechanizing the work through 
the use of SYMBOLIC COMPUTATION. As a result of using this technique, in 
a special-purpose program, the extension of the theory to any order would 
be carried out automatically on the computer. The process flow would be 
as follows: 

I INPUT: Coded differential equation, "order" desired etc... 


Processor: implements the 

algorithm, reorders terms etc. 


Ill OUTPUT: Theory, expressed as a set of formulae, to order 

specified 

Parts I, II, III' will be examined in detail hereunder. At this stage, 
however, the advantages of symbolic manipulation could already be described 
as follows: 

- Automatic development of the theory to any order 

- Mechanization of substitution, transformation, etc. to auto- 
matically condense, simplify, compare formulae 

- Saving of analyst’s time and effort 

- No errors in algebra, given a pre-tested "correct" program and a 


"good" computer. The latter are not minor reservations, of course. 
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but given a sufficient volume of computation, any analyst is 
bound to make errors, in spite of the fact that most of his 
time might have been spent in rechecking algebraic expressions, 


It is apparent in the literature that with the availability of compu- 
ters came a vivid interest in automating the development of analytical 
theories. Considerable effort^ 4 11 to 4 3 4]^ as ^ een invested in build- 
ing special purposes programs and using existing languages for literal 
computations in perturbation theories for the moon, planets and satellites. 

Thorough surveys of existing systems and problems being studied were 

, , ^ . [4-11] , T „ [4-12] 

made by Davis and Jefferys 

Poisson Series, of the form 


<P (x,y) = £[Pj cos(J T y) + Q sin(^ T y)] 


in which 


x is a n-vector of polynomial variables 
y is a m-vector of trigonometric variables 
j is a n-vector of integers 

P. , Q_. are polynomials in the polynomial variables having, 
possibly, negative exponents also. 


have been the object of numerous special purpose programs 


[4-13 to 4-26] 


which aim at economically and efficiently performing the following manipu- 

i .. [4-11] 

lations 

- creation and annihilation of series 

- parsing 


- differentiation and integration 


substitution 
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- numerical evaluation 

Among many contribution, we could mention Barton 1 s^ 4 14 ^ expansion 

of the lunar disturbing function to the tenth order in the LP expansion 

and his attempt^ 4 to develop Delaunay’s canonical transformation 

for the elimination of periodic terms, which ran into severe space problems 

on the computer • This ff space 11 constraint, and the enormous amounts of 

Central Processing Unit (CPU) time needed for list processing are the 

two major difficulties encountered in symbolic computation. Eckert and 
f 4— 1 5 1 

Eckert 1 J used an IBM 620 and a Symbolic Programming System to obtain a 

lunar theory of increased precision. 

Deprit 1 developed an analytical theory of the moon based on Lie 

[ 4 - 32 ] 

transforms , using a set of processors developed for series manipula- 

r 4 - 22 I r4-20l 

tion . Deprit and P^om also used Lie transforms and series 

processors to develop the analytical solution of the main problem in 

Satellite Theory (all gravitational harmonics are zero except J^) . 

[ 4-271 

Carpenter developed a program for automatic computation of general 

planetary perturbations to first order, using Hansen’s theory. Sei- 
delmann 1 J modified Hansen’s method by using an iterative process in- 
stead of a Taylor series. 

To conclude, it can be said that the adoption of techniques or proce- 
dures to economize on time and space results in a restriction of the class 
of problems that can be handled and conflicts with characteristics of porta- 
bility and readability of the programs. The more complex and more specialized 
the system, the more dependent it becomes on specific hardware configurations. 
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4.3 The Choice of Elements 

It is recalled that Lidov's original theory 1 was developed in 

terms of the true anomaly, y, as the independent variable. Yet only 
five elements appear in the differential equations expressing the rates 
of change of the parameters with true anomaly, namely a, e, i, w, ft 
(angles are referred to some plane of reference, not necessarily the 
orbital plane of a disturbing body, but so that ft is defined). An ele- 
ment is missing: it could plausibly be T A (time of perigee at perigee) 
or (mean anomaly at epoch) . The inclination of this element in the 
present theory which thus modifies Lidov's theory, will be shown to be 
critical in t iming the occurrence of perigee in the orbit, as opposed to 
geometrically defining the trajectory. 

In the present study, M^, mean anomaly at epoch, was chosen as the 
sixth element. In Appendix A- 1 , the derivation of the differential 
equation for is given, with the procedure for computing the elapsed 

time from M^. 


Since Lidov, in his work, only used five elements, it implicitly 

amounted to assuming that the unperturbed period of the satellite adequately 

[ 4-2 9 1 

represented the flow of time along the orbit. M. Moe made a similar 

assumption, which is also present in the analog application of her 

r 4 _ 2 .nl , . [ 4-4 to 4- 6 ] , 

work 1 . In Musen s work on long-perturbations , where 

mean anomaly is used as the independent variable, or on short— period 

perturbations^ 4 ” 31 t0 4- 32 ^, in which the eccentric anomaly is used, the 

perturbations of a sixth element were considered. 



4-8 


Roth ^~ 33 t0 ** , in. his attempted extension of Lidov's study, used 

d t 

time t as the sixth variable. The equation for was arrived at by 
extending and truncating the equation relating t to v. From equation 


~ [i + — F cos v - — (1 + £)F sin v] 

dt r 2 ye r ye p t 


or 


^ ^ 2 i 

4^ = -y= [l + — F cos V - — (1 + -)F sin v] 

dv vyp ye e ye p t 


2 2 2 n 

■ 7 = [1 - — F cos v + — (1 + -)F sin v] + 0(ef) 
/yp 1 ye r ye p' t * 


(4.3-1) 


if the forces are of 0(e ft ). Neglecting terms of 0(e A ), the equation for 
the perturbations of time is arrived at by considering the first term in 
the above bracket to be the 2 -body expression (which it is not) , or that 


(— ) _ (^-£) = (—=) [- F cos v + (1 + -)F sin V] (4.3-2) 

dV'pert 'dV 2b /up 2 b Ue r ke P r 


pert 


VS? 


In our opinion, there is a basic flaw in this approach. In the first f j._ve^ 

equations as in Lidov and in our work, Y has been set to unity, whereas in 

the sixth equation, Y has been expanded in order to generate an equation for 

time t. This equation is linear in the forces, and thus allows superposition. 

But nowhere does the feature of double integration appear , which normally 

C ^ * 3 5 1 

accompanies this element in any perturbation theory. To quote Kovalevsky 
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"....We must therefore carry out a double integration 
of the equation giving the semi-major axis before 
being able to obtain the mean anomaly. This is an 
important general result. Whatever method is used, 
no problem of perturbed trajectory can be solved 
in celestial mechanics without carrying out a double 
integration at some stage " 


In the light of this comment, we can recall that _t is not an osculating 
element. Thus in the perturbed motion, the first term in the bracket of 
Equation (4.3-1), which is an "instantaneous analog" to the mean motion 
will be perturbed as compared to its value in the unperturbed case, the per- 
turbation being of order e and noted Oi(e). Thus 

0(1) + OjCe*) + 0 2 (e*) 

where 0 2 (s) corresponds to the two last terms in the bracket of Equation 
(4.3-1). This yields 

At = / 2ir (At)dv 
0 dv 

= [0(1) + 0*(e A ) + 0 I1 (e) ] 
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1*2 

and 0j.(e A ) results from the perturbation of * Thus, to be consistent, 

perturbation of the mean motion should be included, which in turn will 

da 

lead to introducing into the sixth equation a term involving a + dv 
before integrating with respect to v over (0, 2ir). As Equation (4.3-1) 
illustrates, this Roth has failed to include. 

This problem does not arise with osculating element M A (mean anomaly 
at epoch) adopted in the present study. We have, as in Appendix A-l , 



2r 3 . 3 nt da 

l/ 2 r 2 a dv 
yae ^ 



(cos 


. dfi du 
1 d^ + ^ 


) 


It is seen that the unperturb ed mean motion does not appear alone, and that 
the term linear in time accounts for the double integration. 

In Appendix A-l , the sixth quantity used by Roth has been evaluated 
for the J2Q term of the earth's oblateness. The result has been compared 
against a digital integration program (EOLA-TP) run for a high eccentricity 
orbit with perturbations only. Pig. A .1 shows this comparison for the 
time of passage at perigee using either Roth's sixth element or M*, or 
EOLA-TP. Table 4-1 also gives the data from which Fig. 4.1 is plotted. It 
is seen that Roth's elements predict negligible change in the apsidal period, 
which is not true. Using M* gives a close approximation to the digital 
integration. Incidentally, Roth's work does not at all consider the 
sometimes significant effects of oblateness in orbits of high eccentricity. 


4.4 Number of Terms 'Az^ 1 to be Retained 

With the same notation as in Chapter 2, let Az^. designate a change in 
osculating element z over an orbit of the satellite (more precisely , over 
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an interval (0, 2tt) in the true anomaly v). The first subscript, _i, 

is the "order" of the Legendre Polynomial term considered, in the LP 

expansion of the perturbing force F^. The second one, j_, decreased 

by one, the order of the derivative retained, for the term studied, in the 

abed 

Taylor’s series expansion of products such as £1 5 2 ^3 ^ about the posi- 
tion that the disturbing body assumes at the reference time, t ref * 

One remembers that in Chapter 2, expression y 

2 2 _ 1 

y = [1 + — F cos v - — (1 + -)F sin v] 

1 1 ye r ye p t 


such that 

dt _ I 

dv ^ ^ 

never differs from 1 by more than 0*8%, at maximum, in the worst possible 
case, for e = 0,95. It was therefore taken to be 1. Barring other con- 
siderations, the approximation on Y should set a lower bound on the terms 
to be generated. 

Ref. [2-3] as said above, gives estimates of "maximum" amplitudes^ 4 
for the (Az) . To recapitulate the formulae 
a) along the LP expansion, 


< /iL_\ 3±2 2q+5 

~ V ’pk/ q+1 q+3 
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where j is the order of the derivative retained in the term being 
considered. 

On the basis of these estimates, tables of relative "maximum 11 ampli- 
tudes of such terms were drawn up: Table 4-2 lists three such examples, 

for various "large" eccentricities, with the moon as perturbing body. 

Table 4-3 contains an example for the sun. Needless to say, we should 
expect that for the moon, it will be required to much higher values of "i" 
and "j" in the case of orbits of large eccentricities, than might be the 
case for the sun as the perturbing body. An answer as to how many terms 

ought to be retained in a specific case is of great interest. Previous 

f 4- 3 31 

work did not focus on the question: Roth did not set up a precise 

estimate, and Lidov examined the "11", "12" and "21" terms for five ele- 
ments only. Table 4-2 shows that for orbits of e > 0.92, the number of 
terms needed, if the approximation of y were regarded as the criterion 
would be very high (beyond q = 5 in Az^) for the moon as perturbing 
body. Now, despite the use of symbolic computation, the required computer 
time and memory space are very high for high orders of the theory, as will 
be discussed later. In all cases tabulated in Table 4-2, it is seen that, 
for the moon as perturbing body, 

A 51 0.12(A U + A 21 + A 31 + A 41 ) 

Afci £ 0.09CAH + A 2 i + A 31 + A 41 + A 51 ) 

1 being the (normalized) magnitude of the "11" contribution to ^ z tota ^* 

As a tentative cut-off point, the value 0.1 was adopted for the factor 
multiplying the parentheses in the above expression, and the theory 
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developed includes terms 

All> Aj2> Ai3, &21 > ^22 » A31, A4i,and A51 


In the section where the results of the developments are checked against 
high-precision numerical integration, at will be seen that this order of 
the theory is indeed adequate in practice for orbits of e < 0.95. 

To conclude this question, we might remark that the relative effects 
of Sun and Moon, in the "qth' T force term can be listed, as in Table 4-4. 


Let, from TCef . [4-2] , 


F I 5 
q 'max 


\r q 

*f(~) (q+1) 

r d d 


In a perturbation equation, that for — example, this gives 


dft A 
dv p 


p sm 1 


r (q+«(f-) q+2 


sin(w + \>) 


and the relative effects of sun and moon are measured by the ratio 

(Afl) s 

In this estimate, sun and moon have been assumed to describe coplanar 
orbits. The ratio is listed in Table 4-4 for low values of q. The same 
estimate would apply to the other osculating elements. 


y s , M.q+2 
■v — ( — ) H 

U m r c 
tn o 


4.5 The Effect of Earth's Oblateness 

The effects of earth's oblateness on orbits of large eccentricity 

[4-36] 

have been found to be quite significant on natural satellites or 

r 4 - 37 to 4-38] 

artificial satellites, such as IMP-G L • In the present study. 
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-3 

only will be retained, since other harmonics are of order 10 or 
less compared to it* Its secular effect, obtained by integrating the 
corresponding perturbation equation with respect to v, will be superim- 
posed on the effects of the sun and the moon* It is well known that 
there are no secular variations in elements a, e, or i due to J 2 q* 
secular variations in and due to J^q are evaluated in 4.10. 


4.6 Symbolic Integration System 

In the present section, a field of integrals is defined which are to 
be computed in order to obtain closed form expressions for the changes 
of the elements, over one orbital period of the unperturbed orbit. The 
recursive relations involved in the calculation are obtained. An estimate 
Is made of the "explosive 11 growth in the number of terms to be calculated 
for the "Az,. 11 contribution. The results of this study stress the de- 
finite need to resort to symbolic computation for error-free algebra and 
for bypassing the formidable task of hand computation. The elements 
involved in the choice of a particular programming language are discussed. 
The system is described in its various parts. Finally, the relevant 
programs and space and time estimates are given. 

4.6,1 Field of integrals 

It has been seen in Chapter 2 that the integral form ! (I.F.)' to be 


dealt with needs 
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(I.F. ) 5 /*" (At) S S±nUv ■ C ° S - -~ dv 

s,u,v,q def 0 (1+e cos v) «l 

in which 

s,u,v,q are non-negative integers 


( 4 . 6 - 1 ) 


q > 1 ; 

q > u + v 


(4.6-2) 


Integrating (4.6-1) implies a series of recursive operations in 
a finite field of expressions (A 'field' is defined here as a set of 
expressions closed under integration, such as polynomials in several varia- 
bles, polynomials in sine and cosine of an angle, etc.). Basic form 
(4.6-1) will in the process of computation evoke only linear combinations 
of numbers of the set. The element members of the field of expressions 
are not separately tabulated, but the recursion relations hereunder are 
defining each and every one of them. 


4.6.2 Notations 

Let A = 1 + e cosv 


. u v 
sin v cos v 


u,v,q 


H 


. u 
sm v 


u ,q 


= L 

q u,o,q 


v 

cos V 


v,q 


= L 

q o,v,q 
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_ sin v cos v _ L 
v,q " A q l»Vj q 


K 


. u 

- sm v cos v _ ^ 


u,q A q 


u,l,q 


M = — = L 

q A q o,o , q 


= / M x dv = — yy — arc tan [ tan 


and 


2tt 

0 


2tt 

M* 


Let S 5 / H„S . dv 
m 2 m-1 


R x = / M 2 dv 


(4.6-3) 


c = JE -JS. 

- if ^ ~ O 


Vp 3 ' s 3 ' 2 


Note that 


M r 


Sr (At > = 


/ 


yp 


(4.6-4) 


At is the time measured from some fixed reference time, t anc * should 

be considered as one variable, in contrast with the notation ’A' defined 
above, or A in Az^ , designating the change of element z in the "ij" theory. 


4.6.3 A set of vanishing definite integrals (I.F.) 


s,u,v,q 


The reference time, t is chosen as that corresponding to the occurrence 

[4- 1 I 

of apogee along the unperturbed orbit . Hence, At is an odd function 


of v about v=tt , 
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- - < At < - 
2 2 


and 

(4t) v - -(4t) 2j . v 

It is readily apparent that the integrand of (4.6-1) is odd or even 
about it = v, depending on the sum of the exponents of (At) and 
sin v, (s+u). Thus 


r 2x 


(IF) 


s,u,v,q 


(At)' 


. u V 

sin v cos V 
(1+e cos v) q 


0 if (s+u) is odd 
( =#) if (s+u) is even 


(4.6-5) 


and the corresponding sets of (IF)'s can be ignored in what follows. 


4.6t 4 Recursive relations for s = 0 

Let s=0 in (IF) . The case s 4 0 will be treated in 4.6.5. 

s,u,v,q 

In all recursive relations, the goal is to monotonically decrease the 
value of the integer value of q (q - 1) . Since q has to be larger than 
u+v, non-negative integers, u and v will also have to decrease to 0 
or 1. After applying the relevant recursion forms the final results 
obtained is a lengthy primitive expression, evaluated by substitution 
of the limits 2 tt and 0 for V. 


4. 6. 4.1 Integral of L 




Using integration by parts, and the fact that 


sin v , 1 

dv - 


(q-l)e A q-1 


(4.6-6) 
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it is easy to obtain 


L dv 

u,v,q 


1 T _ ■ 2 - 1 - 

(q-l)e u-l,v,q-l (q - l) e 


, L u-2,v+l,q-l 


dv 


+ 


v 

(q-l)e 


L - -dv 
u,v-l,q-l 


(4.6-7) 


with q > 1 . 


4. 6. 4. 2 Integral of K 




0 * O 

By integration by parts, and using (4.6-6) and cos^-v = 1-sin 
we get 


K dv = 


K 


J u,q (q-l)e u-l,q-l (q~l)e 


u - 1 


H u,q-l dV (q-l)e J 


H 


u-2,q- 


with q > 1 . 


4. 6. 4. 3 Integral of H 

u,q 


Using 

cos v - -(A-l) (4.6-9) 

e 

we obtain 


H dv = 


. u-1 , 

sin u-1 


u,q (q-l)e q-1 (q-l)e J u-2,q-l 


K 


.dv 


(4.6-10) 


or 


H dv 
u,q 


(q-l)e H u-l,q-l (q-l)e 2 J 


u-1 


H „-!, q -2 d ’ 


+ 


u-1 

(q-l)e 2 


j H u-2,q-l 


dv 


v , 

l dV (4.6-8) 
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4. 6. 4. 4 


Integral of I 
S 


This computation requires the value of 


I dv 

v-l,q-l 


i .e . 


I _ .dv = 
J v-l,q-l 


f v-1 

COS V 




dv 


v-2 

sm v cos v 

*q-i 


+ (v-2) 


V-3 , 2 

cos v sin^v 


vq-i 


dv 


- (q-l)e 


sin 2 v cos V 2 v 


dv 


Solving for 


I dv, using cos 2 v = l-sin 2 v 

J v,q 


I dv = - 
J v,q 


1 

(q-l)e 


j + ( v -V 

v-2, q-1 (q-l)e 


I . -dv 
v~l,q-l 


+ 


I „ dv 
v-2,q 


(v-2) ' 

(q-l)e j 


1 -3 i dv 

v-3,q-l 


(4.6-12) 


valid for v > 2, q > 2. 


For v = 1 



M 


q-1 



M dv 

q 


(4.6-13) 


4. 6.4.5 


Integral of J 

y»q 



l 

(q-l)e 


1 i + 

v j q-1 


(q-l)e 


J dv 

J v-1, q-1 


(4.6-14) 


with q > 1 
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4. 6. 4. 6 Integral of M 


,M dv = 

J <1 


dv 


( sin 2 v , 

= J—^~ dv 


f COS 2 V 


dv 


sm v 


(q-l)e A q-1 (q-l)e 


COS V 

A q-i 


dv 




fdv 0 

dv + 

.q-2 2 

J 

A’’ 1 J 


fdv 


(4.6-15) 


wherein the identity 

cos 2 v = -^2 (A 2 - 2A + 1) 
e 


has been used. Using (4.6-13) in the expression for 
and rearranging, 

2q-3 


,M dv , e = 1-e 2 , 

J q 


V v = ■ TFUT + 7?i) J • (q-i)qj 


H , dv + 


M ,dv 
q-1 


_q-2 


(q-l)e J 


M 0 dv 
q-2 


(4.6-16) 


valid for q > 1. 

For q=l, we defined 
M dv =S 1 


(4.6-17) 



4-21 


Also, from (4.6-16), written for q = 2 


0 

R = M 0 dv = - - H- - + ~ 
112 e 1,1 e 


(4.6-18) 


4.6.5 Recursive relations for s > 0 . 


The recursive relations for the integration of (IF) 


s,u,v>q 


(At) S L dv - At S L dv - At S 1 M 0 [ L dv]dv 
u,v,q j u,v,q C, J 2 J u,v,q 


From the formulae of the previous section. 


L dv = a Q S, + Fa . (function of type H,I,J,K or L) (4.6-19) 
u, V, q 1 m m 


with the a^s (i = 0,1,2...) being constant. Thus, in (4.6-19), we are 
reduced to integrating known forms, plus a term of the form 


JcAtr^S^dv 


4. 6.5.1 Evaluation of 


S = S -M 0 dv 
m m-1 l 


= S m-1 [ K 2 dv] ~ , M 2 S m-2 [ M 2 dv]dV 


But = M^dv, dR^ = M^dv. Thus 


S R 
m-1 1 


S M^R-dv 
m-2 l 1 


R i r r i 

Vlh - S m-2 2 r + , S m-3 M 2 2f d ' > 
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Finally 

S 

m 


m-i 

Z 

a=l 


(-D 0 " 1 


s 

m-i 



(- I )* 1 ' 1 

(ra-1) ! 


M^R^ m *dv 


(4.6-20) 


Let 


*C = 
m 

Thus , since 


^m-lj 

M i R l dV 


(4.6-21) 


R i = -f H i,i + r 


we get 


t = 

m 




The integrand is expanded by the binomial, formula. For the complete 

evaluation of X , the following relations are also needed: 
m. 


H sj H dv - qfsJ-'Hu H dv]dv 
u,v 1 u,v j 1 1 J U,v 


( 4 . 6 - 22 ) 


and 


q+1 


S? M„ dv = 1 


(4.6-23) 


'1 “1“ q+1 

This completes the set of recursive formulae 
4.6.6 Closure and character of the field expressions. 

The field of expressions is closed under integration under the 
conditions specified in Equation (4.6-2). Ref. 4-2 shows that 


it was necessary to use "mixed" axes for the components of the forces in 

’+■ 

the differential equations of motion, namely (r,t,u); radial, trans- 
verse and normal directions, and (?,^,^): to perigee, normal to perigee 
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in the orbital plane, and normal to the orbital plane. 

It is also evident from the reduction formulae and associated 
definitions, that the field of expressions is essentially not 
polynomial or trigonometric, neither is it "poly-trig”. Transcen- 
dentals such as arctan ^ tan a 2 ) are involved, and the field, 
though finite, is rather general in character. 


4.6.7 Explosion of terms 

The "explosion" of the number of terms under integration is 
briefly reviewed in this section » and given in much more detail 
in Ref. [4-2]. As an example, for the first few orders (q = 1 to 5) 
of the LP expansion, and for each variational equation, we have the 
following number of terms generated, as a minimum : 


q = 1 = number ( 1 ) - 2 

q = 2 ^2 = number ( 2 ) t 15 

q = 3 T 3 = number (3) - 18 

q = 4 = number (4) - 60 

q = 5 T,_ = number (5) - 70 

2 [largest integer in + 2 

following, roughly, a law ^2 

Further, in the Taylor's series expansion of the perturbing 

body motion, the next stage is to multiply each of these terms by 

- 2 and integrate. Each of the above terras again produces at least 
A 

T terras of its own; thus 

q 


T 


total 



2 
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Therefore, both along the LP expansion (q > 1) and the Taylor Series 
expansion (s > 0) , the process "blows up": a large number of terms 
are produced which have to be cataloged and collated during integra- 
tion* As has been pointed earlier, a very significant amount of 
the analytician 1 s time is spent in gathering and checking the results. 
In the approach based on machine symbolic computation, once the pro- 
gram has been checked against known results, the theory can be extended 
to any order automatically, subject to time and memory limitations. 

More important, the recitude of the results is assured, without 
extensive rechecking. 

4.7 The Choice of a Language 

Due to the abovementioned "general" character of the field of 
expressions, we are prevented from using a prepackaged symbolic inte- 
gration program for the generation of the theory. A special system 

[ 4 - 2 1 

for integration within said field was written , in a suitable 

programming language. The latter was chosen on the following con- 
siderations : 

- the language should "match" the character of the field 

- the program at hand uses a large amount of normal, 

- the language should be capable of numerical work 

- it should be generally available, for portability 

Additional factors were: compact representation of elements and 
functions; dy nam i c allocation of storage (space saving; growth of 
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expression not always predictable); recursion, garbage collection 
features; compact storage of lists in a canonical form (to recognize 
common or identical sub-expressions); facilities for rational 
arithmetic (accuracy, recognition of identities). 

Among the languages examined for the task, were ASSEMBLER, 

T 4-'3 9l 

FORTRAN, ALGOL, PL/1, SAC, FORMULA ALGOL, LISP and FORMAC. FORMAC 1 
was retained here for the following reasons: 

- it satisfied most of the above requirements 

- it is a general purpose algebraic manipulation language 

- it has built-in simplification and substitution procedures 
through a canonical form of storage of expressions 

- it is embedded in PL/1, making the arithmetic, control, 
recursion and dynamic features of PL-1 readily available. 

- PL/1 string storage is used for lists not being processes, 
thereby counteracting the storage expenses of FORMAC (double 
word at each node in a list) 

- space allocation can be controlled by an intelligent use of 
a list-erasing instruction 

- the program is readable; it performs algebraic operations 
as easily as arithmetic ones 

- it is widely available at most IBM scientific computer 


installations 
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4.8 Overview of Theory Generator 

Fig. 4-1 shows the main parts of the Theory Generator used to 
generate the extended, modified Lidov's theory contained in what 
follows . 

Equation Generator : it accepts as input the order of the LP expansion 

of the force term (i.e. the highest value taken on by i) and delivers 
as output the variational equations for all six elements considered 
here, with the force components internally generated. 

Preprocessor : this minor routine scans the output of the Equation 

Generator and collects terms which have an identical "operative part. 
The symbolic integration is performed on the equation 

C 1 (At) S f(v) 

with f(v) containing sin v, cos v and A terms, and (not necessarily 
numerical) a coefficient, and preprocessing assures that only the 
"operative” part (At) S f(v) is operated upon only once, with the 
as coefficients. 

Symbolic Integrator : this "core” program accepts as input the varia- 

tional equations, and the order of the Taylor Series expansion. The 
equations (to the order specified) are integrated, and the integrals 
are evaluated between the limits 0 and 2ir in v, true anomaly 

Simplifier : accepts the "new" output of the Symbolic Integrator, pro- 

cesses it through a series of substitutions and simplifications. 
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collects terms and generates, in an internally specified format, the 
final theory. 

The design of the package is such that each program is a "block”, 
written in FORMAC, independent of the other blocks, with communication 
between blocks, or to the user, effected through punched cards or 
files in punched card format stored in mass memory. Printouts of the 
input, and printouts or storage of the outputs at each stage facilitate 
the checking of the flow of information. 

The final output will be the "theory", as a file or a set of 
punched cards. In order to render the theory usable for the user's 
numerical calculations, a small amount of further processing was re- 
quired, more specifically: 

— the replacement of integer fractions by their decimal equi- 
valents ; 

- the definition of "user's variables" to replace common sub- 
expressions appearing in the variation formulae for more 
than one element, which will improve the speed of numerical 
calculation using the theory. 

The formulae obtained on the final result of the theory were 
then incorporated in a PL/1 program, V0LER, described in Section 4.11. 

As regards storage requirements for expressions, they are in- 
herently high, since each operator, operated and associated pointer 
occupies one or two words of storage. Formula manipulation systems 
and user programs are also very cumbersome. Thus, in order to 
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successfully execute any large program, it is of importance that space- 
conserving and -releasing techniques be used whenever possible* In 
FORMAC, this is made possible by the commands 

- SAVE, which stores unneeded expressions on disk; 

- ATOMIZE, which erases a list and releases the corresponding 
space 

Furthermore, the compact string storage feature of PL/1 has been used 
to a considerable extent to economically store lists not in immediate 
use. ATOMIZE has also been used profusely, but SAVE has not been used 
desirably, as it is too expensive in time. Instead, at every stage, 
as soon as an output is generated, it is stored on to a file, or 
punched out, and the space is released by erasing the list. In the 
trade-off between time and working storage, a penalty has often to 
be paid in time, due to the memory space limitations of most digital 
computer systems. This might take the form of integrating again a 
previously integrated "operative part" of which the result had to be 
outputted in a previous run, or of allowing the package to be segmented, 
with the user interacting with the system between blocks. The 
attached risk of error is decreased, however, by visual checks against 
user-or system-created errors. 

4.9 Details on Theory Generator 

Each of the above parts is now briefly analyzed. 

4.9.1 Equation Generator: Program FREQN 

Fig. 4-2 gives a flow chart of the program, and Pef. [4-2] con- 
tains a complete listing- given the recursion formulae. the Legendre 
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polynomials are generated up to one order higher than the required 1 q * 
Their derivatives are computed using the DERIV function of FORMAC, and 
some common symbolic coefficients are generated in the order shown in 
the last box of the flow chart* 

As soon as it is generated 3 each differential equation is expanded* 

d z 

Each rate of variation of an element z with v, — will be a series of 

terms ^ c.f.(v), where as described before the c. are coefficients 
i i 1 1 

and the f^v) are "operative parts". Each separate such term is punched 
out or stored on file. 

Along with the whole FORMAC system and service routines, the pro- 
gram requires: 

- 107 k bytes of memory on an IBM 360 

- for BORDER, 4 (or "q" = 4) , a CPU time of 7.5 minutes (CMU 
IBM 360/67 TSS) 

4.9.2 Preprocessor program COLLECT 

The function of COLLECT, as mentioned above, is to gather terms 
having identical operative parts. This is carried out separately 
for each differential equation, in order to save time by avoiding re- 
peated integration of identical operative parts. A listing is given 
in Ref. [4-2]. 

4.9.3 Symbolic Integrator 

This "block" implements the recursive procedures of Section 4. 6 , 
with the appropriate control, parsing and evaluation routines. A flow 
chart of this block is given in Fig. 4-3. 

I/O ROUTINE: Procedure INPUT 1 

Data : - order "j" of Taylor's series expansion desired 

- output of FREQN, for given q, one term at a time 
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Function: the term is passed to the supervisor for integration 

Output: the final integrated output is printed and punched 

out, or stored on file 
SUPERVISOR: Procedure HIGHORD 
Input : term furnished by I/O routine 

Function: a) examines the terms to be integrated 

b) if it has no operative part, it is returned to 
the I/O as a "constant" 

c) if it has an operative part, the latter is isolated; 
the Taylor index ,j , is updated if necessary; a call 
is issued to the Pattern Recognizer 

d) after completing the integration, it calls on the 
Evaluator 

e) it also monitors integration that might be needed 
during the evaluation process 

Output: result of completely evaluated result is passed back to the 

I/O routine. 

PATTERN RECOGNIZER: Procedure SICODEL, with its internal procedure 

PATTERN . 

Function: The term is examined to see if the relation of Equation (4.6-5) 

is satisfied. If so, the integral is set to zero and sent 
back to the supervisor. If not, the exponents of "sin v' ? , 

"cos v", "A" are determined, and the appropriate integra- 


tion procedure is used. 
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Procedures for the integration: 


r 

- DELDEL: for integrating JM^ dv 

- DELCOS: for integrating 

- DELSIN: for integrating 


J 


H dv 
u,q 


dv = 


L dv 

o »° »q 


L dv 

o,v,q 


L dv 

u,o,q 


- DELSICO: for the integration of the general term 


L dv 
J u > v ,q 


Note that in Pig. 4.3, the diodes indicate direction of calls. 
For instance, DELSIN can call on DELDEL, DELCOS or itself, but not on 


DELSICO. 

EVALUATOR: Procedure EVLUAT 


Function: its primary function is to substitute the limits 0 

and 2ir on v, in the integrated result. Further inte- 
grals, of type S^, which might have to be evaluated, 
are obtained by calling internal routine SSVALU, 

which may issue a call to TTVALU. It also determines 

3. b c d 

the derivative of the perturbing body terms £ ? £ A 

1 2 3 ^ 

The TTVALU procedure integrates terms like 


X 


m 


= [R m_1 M dv 

J 1 1 


in which expression (4.6-18}is substituted for R^. The 
expression is carried out. Control is passed back to 
the Evaluator for substitution of the limits. 
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As can be seen, the Symbolic Integrator is fairly complex. In terms 
of time for integration, it is very much effected by the "explosion" 
of terms to be processed, specially along the Taylor Series Expansion 
dimension, j. As an illustration 

- Space occupied in core, including FORMAC and series routines: 161 K bytes. 

- Time for A^ll^ = 1.5 minutes 

Aj 3 M A = 28 minutes 

A^M^. = 9.5 minutes 

A 22 M a = 40 minutes 

AM, =25 minutes 
31 * 

A,_,M. =4.5 hours 

51 * 

This integration is by far the most time-consuming because of the 
ctcl 

presence of term t — , which necessitates going to 

one order higher in (At) S than specified by the input data n j" (order 
of Taylor Series expansion) „ 

- Time for A]^ 3 3 minutes 

A 13 U) - 11 minutes 

A 22 & = 8 minutes 

=10 minutes 
A 41 e =28 minutes 
A 5 i e = 


1*5 hours 
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The above estimates are made on a TSS environment , and are not 
actual CPU time (the system overhead may be as high as 50% under 
peak load) * 

SIMPLIFIER: Procedure SIMPF6 

Function: to compress the result and print it out in "usable" 

form. The set of simplifications is based on visual 
examination of the integrated output. It mainly 
consists of substitution and removal of common fac- 
tors. The final result is collected according to 
the coefficients dependent on the perturbing bodies. 
Consequently, an extremely compact form of the final 
theory results, as compared to their volume prior 
to the SIMPLIFIER. A listing of SIMPLF6 is given in 
Ref. [4-2] . 

Output: a set of punched cards (or file) containing the theory 

for each element, for each pair (ij). As mentioned 
earlier in Section 4.8, two later subsequent steps 
geared to efficient numerical computations are the 
decimalization of fractions, and the labeling of 
common sub-expressions. 

Space requirements: about 100 K bytes, including FORMAC and 

service routines. 

Time requirements: they are quite significant. For example, 

A 1 3 ^> = 5 minutes 

A 2 2 to - 12 minutes 
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A 3 1 to 

= 2 

minutes 

Ai t ico 

= 2.5 

minutes 

A 5110 

= 3.5 

minutes 


4.9.4 Example 

An example, going from raw data to final output, is given in 
Ref. [4-2] . Times taken at each stage are also tabulated. 

4.10 Results of the Theory: perturbations formulae 

4.10.1 The formulae 

In this section, all six osculating elements a, e, i, (o, are 

tabulated for the following orders "ij" (A_„ z t ^ ie change of 

element z) . 


All A 12 A13 

A 2 i A 2 2 

a 3 1 

A m 


A51 


An important point to mention here is that the formulae were checked 

aeainst those of Lidov^ 4_ ^ , for the first five elements and orders "11", 

[ 4-3 3 ] 

"12" and "21"; and against those obtained by hand by Roth , for the 

first five elements and orders "11", "12", "13", "21", "22" and 31 . (Th< 
sixth quantity he used as a variable is discussed in Section 4.3). The 
agreement was complete, which gives us the highest degree of confidence 
of the recitude of formulae produced by our theory generator. 

For completeness, the secular perturbations of SI, to and due to 
the J£q term are also reproduced. 

Notice that, in the formulae, names have been given to some common 
su b-expressions . They may be repeated and should be strictly associated 
with their definitions within the scope of the "ij" order in which they 


appear at any moment* 
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Formulae A i \z_ 

LP expansion = 1st order ((— ) term) 

k 

Taylor expansion = 1st order ("constant" term) 


Let 



W k a 

6 = — <— ) 




Pi 




B 3 = , 


&4 — ^ 2 ^ 3 ^ * ^5 ” ^ 1 ^ 3 ^ * ^6 ^ * 

e = 1 - e 2 , 

T = time at epoch (perigee), x = period of satellite 
n = mean motion of satellite 


Ana = 0 

1/2 

An e = -15 tt 6 e e £3 


A n i* = ^"Yfi t(5-4e)3 5 cos w - eBt, sin u>], 
e ' 

A n n “ 3 tt — t-tx [(5-4e)B 5 sin w + £$4 cos m] 

e 7 (sin i) 

= 3x 6 e ^ 2 (4Bi - 62 -Be) - (cos i)A^fi 

A U M A = (A u a) - e ^ 2 (A n fi(cos i) + 4 U (d) 

+ Y <5 [ (8 + 12e + 15e 2 )B 6 -(cont'd on next page) 


*The coefficient 6 of A^i is incorrectly printed as (a/p,) 

1 (1 19^ K 

instead of (a/p^) in Lidov's paper' 
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-3(l+12e+22e 2 ) Pi 

-21c3 2 1 


A U T 


-(A n M*)/n 


Formulae Aj 2 z 



LP expansion 83 1st order ( C^) 

term) 


Taylor expansion = 2nd order ((At) 

term) 

Let 





6 = 

V t- p 3 

— (-S-) . n, = mean motion of perturbing body, 

P p k * 

e k 


eccentricity of perturbing body 


V 

&S 

true anomaly of perturbing body in its orbit. 

• ' 

= 

ae k 

— * angular speed of perturbing body in 
d ^i r d «2 E _ 

its orbit 

5 1P 

SS 

d0 k ’ 5 2p d0 k * Sp d0 k 


$ 

“ 

A k ’ *p ■ d£^ ‘ - 34 k e k (sln 9 k> 




2 


01 

£= 

(ZSiSip* + 5 1 *p) 0 k 




2 


02 

S3 

(2CaC2 P + + ^2 < ^ ) p‘ )e k 




3 3 - [(eic 2p + ?i p 5 2 )+ + «iC 2 * p J 0 k 

34 - l<€ 2 C 3 p + ?2p53>t + ?2?3+ p ]0 k 

3s « [(CiCap + CipC 3 )^ + 5i53+ p ie k 

3s - * p $ k 

3 ^6 
Ai2 a = - 2 fi aT [ e ( e+ 4 ) ( 3 2 " J“) 

. i' 

+( 2 e 2 +4e-l)(B 1 ~B 2 )] 

1 B 6 

A j 2 e = “ 35et [^(e+4) (0 2 ~ ~) 

t(f - |)( 3 i- 3 2 )] 

A 12 i ” | 5 T[e 2 + ~ e - cos u - 65 sin to) 

Ai 2 fi = I 6 • , T ; [e 2 + ~ e - |] (B5 cos w+ & 4 sin w) 

8 sxn i 9 3 

~ q 3 b 2 4 

Aj 2 w =-3 6 “(g e + 3 e - e - ^) 3 3 

-Ai2^C cos i) 

a 12 m * - | “ (Aj 2 a ) “ £ ^ 2 (A 12 fi(cos i) + A 12 u) 
+ 5 e 1 ^ 2 t(-y- - 40 e - ~ e 2 )& 3 


A i2 T = -(A 12 M ft )/ n 
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Formulae A i 


LP expansion 


1st order ((”) term) 


Taylor expansion 


3rd order ((At) term) 


Let 


— (— ) 

P Pk 


2 

.. d e, 


V 


dt 


de, 
0 k dt 


angular acceleration of perturbing body 
in its orbit. 

angular speed of perturbing body in its 
orbit 


«i th -th 

■lp “ de. ’ 5 2 p " d 0 * 5 3p de 


dtj> 

de. 


" 3 Vk sin e k 


r pp 


2 

dj> 

2 


= 64 k e^ sln\ - 34^ cos 6 R 



[2»(< p - V + 45^ + V pp 3^ 


+ <2^ lp * + e k 

12*<4 - £> + *52« 2 p*p + Mpp^k 


+ (2C 2 S 2p <t> + t2*J\ 


I2*<5ip5 2p - ClC 2 > + 2* (5i?2p + 5i p 5 2 ) 


.2 


+ £l£ 2 <t> pp ]$ k + [^1^2p + Slp^)^ + e ] i 


[2<f>(5 K ~ K l ) + 2<t> U K + £ £ ) 

1 ™ ' 2 P 3P 23 P 2 3P 2P 3 


.2 


+ 5 e ♦ ]e. + kc c + ? £ )♦ + k e ♦J 0 lt 

2 3 PP k 2 3P 2P 3 2 3 P *“■ 


' 2 *<V 3 P - W + Wtf + « lP V 
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&13 a 

*13 e 

Let C 
D 

Then, 

A 13^ 

* A 13 i 

A 1 3*** 


a x 


1/2 ¥ 


B 3 ( -4e " 


f a 2 + *.* + | e* + f) 


1/2 ' 

e N , 6 t 2 c r 26 _ 5 2 e 2 + 

in < A ” a) + 5TT ! 3 e 4 ,e 


il e* - ¥ e3 - £ e* - 1)B 3 


16 


32 


^5 22 / ? 

<- “ e + (tt 2 


1/2 


|§)e 2 + 2e 3 + | e 4 +~-+ f) 


1 /2 ,4 ,9 ? . Jli. _ 

B 4 e ( 3 e + 32 A 8 ; 


1 S - [C sin to + D cos to] 

2 (sin i) ¥ 


1 — — fi [C cos to - D sin to] 

2 ¥ 


-(A 13 fl)cos i 


+ 1 ^ ^ im- H * + <2e + 3e2 + ! e3 

2 ¥ e 


li) 

3' 


+ B9(-#e + §e-2e 2 -f^e3 + f) 

.2 „3 


+ Be( 4 


TT Z , G 6 6 — u. 1 \ ] 

- - — e + A lt T i / J 


8 3 16 
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o 1/2 

a 13 m * “ 2 a ^ 13a ^ " e (A 1 30 (cos i) + A 13 w) 


+ It Ieez( ' 


39 13 2 

I? - T " - 


3 6 64 ' 


+ M^T + e <l * 2 ~ 13 > 


+ e 


2 <§ " 2 


64 


• ,5 2 39 . , 143 9 2. 


‘8 " 16 v 3 2 


143 9 _2 , 

3 


+ e 2 ( i29 _ 17 * 
+ 6 ( 64 4 * } 


13eB-f eM] 


A13T = -(Ai 3 M a )/h 


Formulae A? \ Z_ 

r 2 

LP expansion = 2nd order ((~) term) 

Taylor expansion = 1st order ( "constant” term) 


a 4 

6 = ~ (~ ) , 

■f 

♦ = A. 

v p k 

k 

«i = C i<}> , 

“2 = ^2^ » a 3 = £ 3 $ > 

3 

3 2 

Yi D £ 1 $ » 

Y 2 “ £ 2 ^ » Y 3 = ^ 1 ^24> » 

2 

2 2 

Yi* = ClS3<J> , 

Y5 = ZlZs'b* Y6 = ^2 C 1 ‘f’ » 


Y 7 d C 1^2^ 3 ^ *- 
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A 21 a 

A 2 ie 


= “ ir<Se [ (7-6 s)y 3 + e Y 2 

o 


- (7 - 3e)-f] 


Letting, for brevity, 


a 3 


[(7-4e)y 4 + cy 5 - (7-30“ ] 


ley. 


we write, 
A 2 1^ ~ 


- f * 5 T75- 7STY) l c sln “ + D cos 

e 


A 2l i " - ^ ir <5 (C cos w - L sin w] 


1/2 


A 2 l * 11 


13-V6-Z— [<! - c)Y, * (3 - r)Y 6 - —(13-90)-i 2 l«(=o s » 

2 e 4 _£ 24 20 


a 21 M * 


2 21 (A 21 a) - e l/ 2 (A 21 fi(cos i) + A 21 o>) 

Z 3r 


+ — tt 6 [45e£Y6 +(8+21e+24e 2 +52e 3 )Y 1 
8 


- | a . (8+36e+24e 2 +37e 3 ) ] 
5 1 


A 21 T = - (A 21 M )/n 


"'The term underlined is misprinted as Y 14 in Lidov’s paper 
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LF expansion 


Taylor expansion 


2 nd order (( 7 -) term ) 
k 

2nd order ((At) term) 


Let 

6 

5 

IP 


t <f-> 4 ■ * - \ > » P ■ - 44 k e k (sin V 

v H k 

«, _ «a .111 

d0 k ’ ’ C 2P = de k * C 3P d0 k 



iik 

dt 


angular speed of perturbing body in its orbit. 


°1 ^iP* + 5 iV® k 

“2 " + 5 2 *p )e k 

«3 - <V + SV> 

Yi - '<35^ ip + + «iV® k 


y 2 = 

ok 2 k. 

2 - 

+ 

2P 

5 3 

y\ 



y 3 ■ 

(2c 5 

i 

lp*2* 

+ 

5 ? 5 ZP* 

+ 

Z 2 t. 4 )9. 

s l s 2 p k 

Y4 = 

(25 5 

1 

ip £ 3* 

+ 

E 2 £ 6+ 5 

** 1** 3p v 


Y5 " 

(25 5 
2 

K * 

2P 3 

+ 

4 

2 3P 

4* 

«!«.V 5 k 

Y 6 " 

(25 5 
2 

x ♦ 

2P 3 

+ 

5*5 « + 
2 IP 

+ 

«2 5 lV 5 V- 

Y? = 

[(5 
. IP 

5 C • 

. 2 3 

+ J 

'.Vs 

+ 

VsV* 


+ Cl^2^34* p ]0 k 
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A 2 2 a 
A 2 2 ® 


Letting 

C 

D 

we have 
A 22 $l = 

A 2 2 ^ “ 


2t 6 a[—£ eeY 6 — ^p(12 - 6e + 36e 2 + “ e 3 ) 


+ 8 Y 1 ^ “ 1 e + 12e2 + 6e3 ^ 


, E . 615 2 645 v 

tS ; !< 3Cte + '5T e " 6T )eY 6 


, y 3 ,3 2 _L 3,3. 

+ a^- ge + ^e +^-e +-) 


+ Y ( 155 25 2 _ 15 3 _ 5., 

+ Y l ( 64 e 2 6 4 6 2 )] 


2ae 


(A 22 a) 


e " 25eZ “ e3 " 5 


r r 3 ,3 0 j. 21 3 * 3 . 

[a 3 (- - e + - e 2 + — e 3 + -) 


64 


x. z 255 

% < 6 T e - 2 


25 e 3 - ^ e 3 _ |) 


. , c 135 . , 

*^ 5 ( " 5 ~ 64 7 e)e] 


t 6 . 


(sin i) 


[C sin a) + D cos w] 


T3[C cos w - D sin to] 
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A 2 2“ 


3 , 3 


111 


'l M-! + ! e+ 64 


^ 9s 21 4 \ 

j 2 - l e 3 - -rr e 4 ) 


16 


+ v|?-¥ • - S 1 - 1 + ? •’ * «•'> 


+ 77 Y, (~9 + 6^e + 36e 2 )e] 
64 2 


- A 22 ft(cos *■) 

a 22 M * " | |(A 22 a) “ e 1>,2 (A 2 20(cos i) + A 22 <o) 

1/2 nS 

+ 6te 7 [y 2 E(20 + ^ e) 

+ ^ y 3 <16 - y 1 e + 80e 2 + 24e 3 ) 
+ 9 a 2 (-2 + | - 2e 2 - e 3 )3 


A 22 T « - A 22 M*/n 


Formulae A^ jz_ 

r 3 

LP expansion = 3 rc ^ order (( ) term) 

k 

Taylor expansion = 1st order ("constant" term) 


Let 
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2 

» £i4> » 

2 

a 2 D £2$ * 

a 3 “ 

«4 “ ^2^3^ > 

a 5 15 » 

a 6 = <}> 

4 

$1 ® £i4> > 

3 

82 D £l£2^ » 

2 2 

83 “-5lS2$ 

3 

$4 “ C 1^2^ » 

4 

85 = £2^ » 


3 

Y 1 B SlM 1 

2 

Y 2 “ £l?2^3^» 


2 

Y 3 = ^1^3^* 

3 

Y4 “ £ 2 S34> » 



A 31 a 


*31 e 


Letting 

.. . -C 


0 

■i2W /2 e[303C2 + e 2 > 

o 

- 73z C2e 2 + 1) - 7e8 4 3 


[(630e 2 + 105)- 


8 


,615 
^ 8 


2,135 4 

e + — r- e 4 


+ 


15 

2 



+ (- 


315 


e 2 + 105 


+ J f>T> 


D 


Then, 

A 31 0 


Y 

e 2 + ^)a 4 + (630e 2 + 105)-j- 

. 105 , 

+ -g- cyj 

1/2 

6tt v - -r r [C sin u> + D cos u] 

(sin i) 


A 3li 


Sire 


1/2 


[C cos ui - D sin w] 
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A 31 w = - A 31 C2(cos i) 

+ ufie 172 ^ e & 3 - ^ e 2 6 3 - (540e 2 + 615 )-j 

+ (105e 2 + 

+ ■^•(6e 2 + l)ct2 “ “ 

.45 j ± 15> , 

+ (“g e 2 + ” )a 6 ] 

A 31 M* = " (A 3 ia) -e 1/2 (A 3l fi(cos i) + A 31 w) 

+ 6 it [- ”■ 3 j(| + 8 e + 21 e 2 + 8 e 3 + 20 e 4 ) 

- | a 6 (8 + 12e + 37e 2 + 12e 3 + e 4 ) 

. 11 | 5 . £ 2 65 + 1|5 ct 2 (4 + 3e 2 )e + g~ ai{lO + 48e(l+e 2 ) 

+ 3^1 e z + 87e 4 }- 8 3 (1 + 5c 2 - 6 e 4 )l 

A 31 T = - (A 3i M*)/n 

r 4 

4 th order (( — ) term) 
k 

1 st order ("constant" term) 


Formulae A 4 \Z_ 

L? expansion “ 

Taylor expansion = 


Let 


6 

♦ - \ 
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For economy of notations, only the following symbols will be 


explicitly defined: 



2 

4 

6 

<*1 “ £1 » 

a 2 - £1 

. ct 3 = £j » 

2 

4 

6 

8l 1=5 C 2 » 

3 2 - C 2 

» 83 = £ 2 > 

y 1 0 ^iC 2 . 




Note that these symbols are not identical to, or consistent with 
those used for Ajj, A 21 , A 3 1 ; that <j> appears as a factor in the 

products and that £i» £ 2 , £ 3 * also appear explicitly. 

Ai,ja — 0 

A 4 ^e = m 6 e ^ 2 <j> [- 0 ^ 2(2 + 23e 2 + 8 e 4 ) 

+ a 2 £ 2 (l + I 6 e 2 + I 6 e 4 ) 


+ -g— £ 2 (8 + 20e 2 + 5e 


+ aiSi£ 2 (l + 8e 2 ) 


“ ^ 8lC 2 (2 + e 2 )e 


+ 1 |f e2 M 2 ] 


4 ) 


Letting 

C = - “■ (8 + 20e 2 + 5e 4 )~ 

+ 2205(2 + 7e 2 + 2e 4 )^ a x 
- 6615(1 + 2e 2 )ea 1 8 1 


(cont*d. ) 



4-49 


<5 + 20e 2 + 8e 4 )| a 2 


2205 


ec8 2 + 735(2 + e 2 )e8i] 


[735(2 + e 2 )eY! - 2205e£Yi6i - 2205(1 + 2e 2 )eYia 1 ) 


we have 


A 41 ft = 6 ti 


}/ 2 


(sin i) 


£ 3 <t>[C sin co + D cos to] 


1 /2 

A 41 i = 6ire ' ^ 3^ [C cos 00 - D sin co] 


A4 


§§£ e l7 % [-735 (-2 + 3e 2 + 


+ (-1 + 5e 2 )£?!e 2 


- ^ (8 + 60e 2 + 25e 4 )Ci 


- 2205(1 + 3e 2 - lOe 4 )?^^ 
+ 735 (2 + 21e 2 + 10e 4 KiC»i 


2205 


(1 + 12e 2 + 8e 4 )Cict 2 ] " (cos i)A 41 « 


& 41 m * ” |a ^ Al +l a ^ “ e ^ 2 ( A 41 fi (cos i) + A 41 co) 

+ iT6(j>[- ^“~(2 + e 2 )eeCi3i 


+ H (64 + e + 320e 2 - ~ e 3 ) e Ci6 2 


335 


95 


(cont'd.) 
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+ -J(12 + 122e + 120e 2 + 335e 3 + 6Qe u + ~r- e 5 )Si 
16 H 


- 63(3 - + 15e 2 -20e 3 )eCicti8 1 

- (16 + 102e + 160e 2 + 477e 3 + 80e 4 + 166e 5 )5iai 


+ (43e + 96e 2 + 268e 3 + 48e 4 + 136e 5 + T^)Cia 2 

64 D 


. / 3 4 7 e 3 4. ,945 _ „ zpj:> r „ fi 

+ (- 5 e ~ 77 e " T ~ o) (“5~ ^ 1^2 ~ L Sl a l&l' e 


2835 


8 " 15 8 3 2 


+ (| e + ±f) 5l B 2 - ^ 5ia,6i)e-] 


2835 


8 


15 2 


Ai+lT = - (A4iM A )/n 


Formulae for Ac; i z. - 

r 5 

LP expansion = 5th order ((~) term) 

r k 


Taylor expansion = 1st order (”constant' , term) 


With the same notation as for A 41 , let 



- 

2 

4 

6 

<*1 » 

Cl i 

“2 ** Cl > 

a 3 ” Si 


2 

4 

6 

6 1 = 

> 

&2 " ^2 * 

63 *= C 2 


Yl “ £ 1 * 



^ 51 a 

A 51 e 


3 /2 if 31185 


«5tree <J>[ - 


64 


E Yl &2 


15 


- 315(| I> 


2 )I> 

V c 


+ 945<£ e 2 + \ ^ + ?)^i 


32 


4 e 


c 3 i 15 

- 2079(f e 2 + | e 4 + 


64' e 


+ 2835 <^ + jWjBj 


- 10395(|- + -|j)v 1 iiiB 1 ] + 


2ae 


(^ 5 ia) 


Letting 

C 


[3|65 (e2 + i )£5i62 


+ if 1 (I 1 e 2 + e" + 5e 6 + 1)— 


9 A 5 / 23 o ^ o 

- — (g- e 2 + e 4 + I )Ci3i 


4 'V 


(45e 2 + 80e 4 + 16e 6 + 2)£j 


+ -~^(2e 2 + 2e 4 + |)£ia 1 e l 

+ —• (120e 2 + 240e 4 + 64e 6 + 5)K\— ] 
64 £ 


[“ 3^ ( 2 + e 2 )eC2Si + t 2 h e 
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+ — | <20e 2 + 5e 4 + 8H 2 


3465 , 2 . I*. 


+ ^ (e 2 + | )eS2ai3l 


■ 

945 ,23 2 j. 4 , 1\ f 

- — (g“ e z + eV 5 )C 2 a l 


+ (2e 2 + 2e 4 + |)£ 2 a 2 ] 

we have 

1/2 

A5 jfi = 

6ir '^' s ^ n jj £3$[C sin w + D cos 

Asj! - 

Site ' C 34* [C cos w - D sin tu] 

A51&J - 

c 1/2 6ir^[-°^y-(l - 4e 2 )a 1 8 2 e 


2835 / o , i \ n 
+ -gj- (e 2 + l) 82 e 


+ ■ (10e 2 + 8e 4 - 7)016! 


+ . (75e 2 + 20.e 4 + 26)^ 

J * 


+ ^gf^C-iee 4 + 5)a 2 8i 


- “• (160e 2 + 48e 4 + 45)a 2 


693 

8 


+ (20e 2 + 8e 4 + 5)a 3 


- (10e 2 + 5e 4 - 4)6i 

+ 15 35 2 _ 35 4 „ 7) 

^ 8^2 8 J 


3465 2 

64 


e 2 3 3 l 


- (cos i) 
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A 51 M a = \ \ ( A 51 a ) -e (A51Q (cos i) + A 5 jio) 


+ [- 


1575 

128 


(20e 2 + 5e 4 + 8)cg 1 


14175 

64 

51975 

128 


(23e 2 + 8e 4 + 2)ea 1 3 1 
(16e 2 + 16e 4 + l)ea 2 Pi 


, 14175 /n , % 9 r 17325 

+ ~i2s~ (2 + e )e 32 ■ ~nr e 33 


+ & e + e 2 + ^ e 3 + ^ e 4 + £ e 3 


64 


2385 6 15 

128 2 ; 


- (24e + 135e 2 + 80e 3 + —■ e 4 


+ 24e 5 + 46e 6 + 6)04 

+ (3e + — ^ — e 2 + 10e 3 + 30e 4 + 3e 5 + 7e 3 + -7) a 2 

4 32 lb 

- ■£?§• (12e + 45e 2 + 40e 3 + 120e 4 + 12e 5 + 38e 6 + |)a 3 

16 o 

- (8e 2 + l)e 2 ai 62 l 

A 51 T = - (A 51 M A )/n 


Formulae Az due to Oblateness 

The only term in the oblateness potential considered in the 
following formulae is 
Let K = J 2 (^) 2 

where R = equatorial radius of earth 

p = parameter of satellite orbit 


ae 
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= - 3Ktt (cos i) 

= 6Ktt( 1 - ^ sin 2 i) 

= 3K - (1 + e) 3 (1 - 3 sin 2 i sin 2 m) 


Aft 
Aw 

AM* = 

AT = - (AM* ) /n 


4.10.2 Auxiliary formulae relating to the perturbing body 

It was explained earlier that all terms which depend on the position 

of perturbing body "d" are calculated at a time, t corresponding to 

the occurrence of apogee along the unperturbed orbit. 

Keplerian orbits are adopted as models for the sun and the moon. At 

any time, the orbital elements are calculated using the mean elements at 

epoch 1900 Jan 0.5 and their secular variation 1 . The mean anomaly of 

the perturbing body is similarly calculated. 9^, which is here the true 

anomaly of the perturbing body in its orbit, is computed using the formulae 

of elliptic expansions 1 . 9°, angular velocity, and 9^, angular accel- 

k 

eration, respectively , are calculated by taking the time derivatives of 9^ 

in terms of the mean anomaly, the mean motion of the perturbing body being 

known. The coordinates of the perturbing body in its orbit are also 

[4-4i 1 r k r k 

calculated using relevant formulae for — cos 9, and — sin 0, . 

a k R a k R 

Let [T ] be the transformation matrix from the (P. , 0. , 1 ) system 

r k k 

to the (P, Q, 1 ) system of the satellite. [T ] is independent of 0, . If 
H TT K 

£l> S 2 » 5 3 are the derived director cosines of the unit to the perturbing 
body in the (P, Q, l n ) system, then 


f N 


r n ‘'l 

COS 0. 

k 

c 2 

= [T f ] 

sin 0, 
k 

.53; 


. 0 J 
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dZ. 

1 

and if. as in 4.10,1, ~ -tt— (i = 1,2,3), 

’ ip d6 


f > 

Cl 


-sin 0, 
k 

C2 

= [T ] 
r 

cos 0 t 
k 

.C 3J 


. 0 j 


Similarly, with £. = [?.], (i=l,2,3), 

X PP dfc) k iP 


5 1PP ' 5 1 

C 2pp = ” C 2 


C 3 PP " ? 3 


In other words; second and higher order derivatives of § 1 , Kz* with 

respect to 0^ can be written in terms of ^ 1 ,^ 2 ? £3 or £ip* ^2pj ?3P* 

4.10.3 Some comments 

It is readily apparent that at higher orders of the LP expansion, the 
formulae became increasingly longer and more complex. The formulae might be 
condensed a little by recognizing common subexpressions in a hand translation 
of the formulae. The chances of error, however , might in final analysis far 
outweigh the improvement in computer time. 

With regard to the perturbations in oblateness, note that they depend 


on factor K : 
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A w. 


An 


2 

^ K=J 2 * 
P 


1 

2 

e 


and 


K 1 

AM . ^ — 'v — 5- 
* C e 3 


In orbits of large eccentricity ,e is of the order of 0.1, and changes quite 
significantly overtime. Thus oblateness perturbations are quite sensitive 
to inaccuracies; the next section will illustrate this problem of small 
divisor. 

Finally, as long as any one plane, say the equator, has been used as 
the same reference throughout, the perturbations on the angles, as well as 
those on a and e , which are due to the sun, the moon or the oblateness, 
respectively, can be summed up directly. 
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4.11 Numerical Verification of the Theory 

The "closed form" theory describing the changes in the orbital 
elements, z,, to orders "11", "12", "13", "21", "22", which has been 
obtained by non-numeric computation, is expressed by the formulae 
of Section 4.10 and is implemented in program V0LER, was checked by 
comparing its predictions to those of high accuracy numerical inte- 
gration programs (NASA r s ITEM and C-MU EOLA-T). The present section 
discusses those verifications and analyses the significant gains 
realized in computer time and the level of accuracy achieved. 

4.11.1 Program V0LER 

A program called V0LER (for evolution of Orbital eLements in 
high Eccentricity oRbits) has been written in PL/1, which uses the 
theory of Section 4.10 to predict the evolution of orbits of satellites 
perturbed by the gravitational effects of the sun, the moon and the 
earths oblateness. A flow chart follows, which lists the names of 
all procedures in the program. A brief description of these is given 
hereunder . 

VOLER . 

Procedure VOLER is the main calling procedure. It initializes 
structures and arrays, reads in data and calls all the major proce- 
dures. Input and output are controlled by this procedure. 

SUN . 

Procedure SUN computes the position of the sun at any given time. 
The model used is based on the mean elements of the sun at epoch 1900 
Jan. 0.5 and their secular variation. The procedure calculates the 
position of the sun in the equatorial system of the earth; and then 
computes the array £ of the projections of the unit vector to the sun 
on the orbital axes of the satellite. The array £p, being the deriva- 
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Start 

4 

Read in Data 
4 

Compute Initial 

Orbital Elements 
4 

STABILITY 

4 

Compute Oblateness 
Perturbations 
4 

Compute Lunar 
Coefficients 
4 

Compute Lunar 
Perturbations 
4 

Compute Solar 
Coefficients 
4 

Compute Solar 
Perturbations 
4 

Update Orbital 
Elements 


Procedures 


no 


VOLER 

SUN 

MOON 

TALON 

COORD 

ORBIT 

ECLEQ 

EQORB 

MOONEQ 

CALNDR 

REDUCE 

PERTRBN 

SUMMER 

OBLATE 

UPDATE 



Flow Chart of Program VOLER 
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tives of the array £ with respect to the true anomaly of the sun in 
its orbit, is also computed. 

MOON . 

Procedure MOON performs the same functions as procedure 
SUN but for the moon. The arrays £ and £p, for the moon, are also 
computed. The model for the moon is again based on epoch data at 
1900 Jan. 0.5 and the mean secular variation of the elements. 

TALON . 

Procedure TALON computes the longitude of nodes of the satellite 
orbit , given the geocentric , equatorial latitude and longitude of the 
perigee, the direction of launch (north or south), the time, and the 
sidereal time at Greenwich at 0.0 Hrs. U.T. on Jan. 1 of the year of 
launch. The direction of launch, along with the latitude indicates 
whether the satellite is approaching or leaving the ascending node, i.e., 
the intersection of the positive nodal line with the orbital plane. 

It is always assumed that the satellite is injected into orbit at 
perigee (or that reduction to perigee has been effected elsewhere) . 

COORD . 

Procedure COORD is called by procedures SUN and MOON. It computes 
the coordinates of the perturbing body in its orbit, given its mean 
anomaly and eccentricity, using formulae of elliptic expansions. 

ORBIT . 

Procedure ORBIT computes the components of the P, ^ and 1 


axes 



of the orbit of the satellite in the geocentric equatorial system. 

> 

ECLEQ, EQORB , MOONEQ . 

Procedures ECLEQ, EQORB and MOONEQ transform co-ordinates from 
one reference system to another. ECLEQ transforms from ecliptic to 
equatorial; EQORB from equatorial to satellite orbital axes; MOONEQ 
from the moon's plane to the equatorial system. 

CALNDR. 

Procedure CALNDR converts the time, day and month of launch to 
an equivalent number of days since the beginning of the year of 
launch. This procedure is also called when a satellite orbit decays; 
it then calculates the day, month and year of the collapse of the 
satellite given its lifetime. 

REDUCE . 

Procedure REDUCE normalizes angles to a value between 0° and 

360°. 

PERTRBN . 

Procedure PERTRBN is the major procedure which contains all the 
theoretical results, obtained as formulae and presented in Chapter A. 
It calculates the perturbations of the osculating elements of the 
orbit over one orbital period. 

SUMMER. 

Procedure SUMMER is called by PERTRBN and merely sums up all the 
perturbations of various orders. 
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OBLATE . 

Procedure OBLATE computes the secular oblateness perturbations 
due to the J^q term on the satellite orbit, using the formulae in 
Chapter 4. 

UPDATE . 

Procedure UPDATE updates the orbital elements by adding the 
total perturbations to their initial values. 

The time requirements are: about 0.6 sec of CPU time (IBM 

360/67, TSS , version 8.1) to compute over one orbital period. Thus 
the computing time is inversely proportional to the orbital period. 
As an example 


Orbital Period 

CPU time/j 

rear of orbit 

2.5 days 

90 

sec 

4.0 days 

55 

sec 

5.0 days 

43 

sec 

6.0 days 

37 

sec 


The above range corresponds to 0.92 < e < 0.95, with a low initial 
perigee (h^ ^ 200 km.). 

This should be compared to the time taken for the digital inte- 
gration of the equations of motion over one orbital period by ITEM, 
EOLA-T... Typically, for the higher eccentricities, the CPU time 
might be of the order of 10 to 20 minutes per year of orbital evolu- 
tion. Therefore, in a rough sense (since comparisons ought to be made 



on the same computer, using the same input/output procedures etc.), the 
savings factor of VOLER is in the range of 10 to 60, the high factor 
applying to larger eccentricity cases. This, obviously, is obtained at 
a cost in accuracy, but this cost is often perfectly tolerable for 
many purposes in mission anslysis. 

In the following, a few significant examples are described. For 
a more exhaustive treatment, the reader should refer to [4-2 ]• 


4.11.2 Type of orbits studied in the examples 

Two major parameters characterize the examples studied 
1) Eccentricity we consider "large" eccentricities defined here by 
0.9 $ e i 0.955 

Note that, in the "11" Lidov's theory, it is shown that 


'll 




Pk 


and if we compute the ratios 

(Ae) e=0.945 / Ae) 0.92 = 2,9 


Therefore, it is seen that the perturbation increases by a 
factor of about 3, even over the "small" range considered. 

2) Inclination: both planetary-type orbits, i.e. having small in- 

clination on the orbits of the perturbing body, and orbits quasi- 
normal to the ecliptic (such as for IMP-G) will be considered, 
with data considered by NASA for specific satellites, and in one 
case post-flight data of an actual satellite. 



2) EOLA-T. It is a numerical integration program, developed at 
C-MU under this grant, using a method of variation of para- 
meters (conventional osculating elements), with time as the 
independent variable. Its detailed features are given in 
Chapter 6. Again, for the purposes of the comparison, EOLA-T 
was run only with the sun, the moon and oblateness (the latter 
only where indicated) . 

4.11.4 Some examples treated 

A. High inclination orbit (IMP-G type orbit) 

These are described best by the tables below comparing the re- 
sults of VOLER with those of a high accuracy numerical integra- 
tion, which can be characterized as follows 
1A: "medium high" eccentricity ^ 0.93 with OBLATENESS 

IB: "medium high" eccentricity ^ 0.93 without OBLATENESS 

The agreement appears very good, with errors of the order of 
2 days/year (0.55%) on the timing of perigee, 0.04%/year on 
the eccentricity. 



TABLE 



Tni t inl Orbital Data: II'F-G - 

Example 1 



1A 

IB 

* 

H-Per : 


405.62 

402.78 

a : 

. 

95,804.57 

94,940.95 

e : 


0.929191 

0.928577 

i : 


86,8665 

86.8659 



105.8008 

105.8045 

R : 






-160.0022 

-159.9953 

w : 




t. 

Year 

1969 

1969 

in 3 





Month, Day 

June 24 

June 24 


Hour (UT) 

17.96431 

17.96448 


7% 

In all tables in this chapter, the following abbreviations 
are used: 

N : 

Orbit number 

t : 

Time since injection (days) 

R-Per : 

Distance to perigee (Km) 

H-Per: 

Height of perigee (Km) 

.a : 

Semi-major axis (Km) 

e : 

Eccentricity 

i : 

Inclination (deg) 

R : 

Longitude of nodes (deg) 

u : 

Argument of perigee (deg) 

t . . : 

Time at injection into orbit at perigee 



TABLE 


Comparison of Numerical Results (IHP-G - Example 1A ) 


* 



Theory 

N. I. 

Theory 

N.I 

N 

53 

53 

107 

107 

t 

179.28 

178.69 

362.70 

360.77 

R-Per 

8,128 

8,123 

9,447 

9,430 

H-Per 

1,750 

1,745 

3,069 

3,052 

a 

95,746 

95,412 

95,690 

95,132 

e 

0.91511 

0.91486 

0.90128 

0.90087 

i 

86.55 

86.41 

86.51 

86.46 

a 

105.23 

105.11 

104.87 

104.83 

W 

200.04 

200.04 

201.57 

201.47 

No lifetime 

figures available 




ft 

N.I.: Numerical Integration 

(Program ITEM) 




TABLE 

Comparison of Numerical Results (TMP-G 

- Example IB) 


— 

Theory 

* 

N. I. 

Theory 

N.I. 

N 

53 

53 

107 

107 

t 

178.66 

178.69 

360.66 

360.78 

R-Per 

7,745 

7,763 

7,970 

7,968 

H-Per 

1,367 

1,385 

1,592 

1,590 

a 

94,908 

94,927 

94,831 

94,844 

e 

0.91840 

0.91822 

0.91596 

0.91599 

i 

86.47 

86.46 

86.90 

86.78 

n 

105.79 

105.78 

106.14 

106.06 

u 

202.98 

203.05 

206.42 

206.59 

* 

• N.I. 

: Numerical Integration 

i (Program ITEM) 
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B. Low Inclination orbits 

The orbital inclination on the orbital planes of the perturbing 
bodies is relatively low . Considering, as one example among 
many listed in Ref. [4-2], a "very high" eccentricity orbit 
such as that of IMP-I (e„ - 0^4 3). 

The comparison between the results of this theory, through 

VOLER, and of a numerical integration are given graphically 

(Fig .4-4 to 4-7 ) and also in the table hereunder. It is seen 

that an accumulating error is present , which however does not 

grow to be very large at the end of one year: about 1% in 

the timing of perigee and in the distance of perigee due to the in— 

2 

accuracy in a and 1-e . Yet, the errors are not unduly 
large, and the overall trend is sufficiently well captured at 
a savings in computer time of the order of about 50. 
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TABL E 


Initial Orbital Data: IHT-I - Example 


H-Per 

236.28 

a 

115,067.60 

e 

0.9425169 

i. 

28.7763 

ft 

216.0352 

U) 

302.3777 


t . . 

inj Year :1971 

Month, Day. ....... .March 13 

Hour (UT) 16.00 
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TABLE _ 

Comparison of Numerical Results (IMP- I - Example ) 



Theory 

A 

N.I. 

.Theory 

N.I. 

N 

40 

40 

80 

80 

t 

179.24 

177.83 

359.13 

355.7 

R-Per 

14,535 

14,256 

22,858 

23,116 

H-Per 

8,157 

7,878 

16,480 

16,738 

a 

115,103 

114,186 

115,047 

114,240 

e 

0.87372 

0.87515 

0.80131 

0.79765 

i 

37.54 

38'. 81 

44.34 

43.36 

n 

193.18 '• 

193.13 

184.78 

186.48 

a) 

324.4 

324.38 

334.53 

332.70 

* 

N* !• : 

Numerical Integration (Program ITEM) 
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4.12 Conclusions 

The present chapter demonstrates the power of non-numeric computa- 
tion to generate a closed-form theory (from perigee to perigee) for 
high eccentricity orbits. An extended, modified Lidov's theory has 
been developed and implemented in a numerical program VOLER, which 
simply evaluates the values taken by the symbolic expressions obtained 
after one satellite revolution. This method seems to be ideally suited 
for calculations in a mission analysis, where requirements for ex- 
tremely high accuracy might be treated for the low computer time and 
ease of use of the present approach. 
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TABLE 4 . 1 

Comparison of Roth's Sixth "Element", Our Sixth 
Element and Numerical Integration 

Initial Orbital Period of Satellite = 5. 04680 days 
Orbit No. Time at Perigee (in days) by 



Roth's Sixth 
■ "Element" 

Our Sixth 
"Element" 

EOLAT 

1 

5. 0468 

4.9621 

4.9629 

5 

25.2339 

24.8073 

24.8148 

10 

50.4677 

49.6068 

49.6294 

15 

75.7015 

74.3987 

74.4442 

20 

100.9354 

99.1831 

99.2599 

30 

151.4029 

148.7276 

148.8918 

40 

201.8705 

198.2434 

198.5261 

50 

252.3382 

247.7310 

248.1615 


6 1 
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TABLE 4-1 

Orders of Magnitude of Terms: LP and Taylor Expansion 

v "Extremely high" eccentricity 

e = 0. 95, a/p^ = 0. 342 for the moon, (rn/n^) = 0. 588 

for the moon 


j = 

1 

2 

3 

q = 1 

1.0* 

0.6 

0.17 

2 

0.9 

0.54 

0.15 

3 

0.74 

0.44 

0.12 

4 

0.58 

0.35 


5 

0.44 




6 


0.32 
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TABLE 

Orders of Magnitudes of Terms: For the Sun 

-3 

0.95, a/p, =0.88x10 for the sun, (7tn)/n =0.05 

K iC 

A . 

qj 

j= 1 2 

i = l 1.0 0.05 

2 2.3 x 10~ 3 


TABLE 4-4 


Ratio of Perturbations Due to the Sun and the Moon 


u /u. 
sun moon 


hh , , c 

= — c = 3 x 10 x 81 


P P 


m 


r /r 
moon sun 


= 2.43 x 10 

= 3.844 x 10 5 /(1.5 x IQ 8 ) 


= 2. 56 x 10 


Let R = ratio of the effects of sun and moon for the 
^ q-th LP force component 

Then, = 2. 43 x 10 ? x (2. 56) 3 x 10 _9 


=* 0. 39 

R 2 = 2. 43 « 10 7 x (2. 56) 4 x 10' 12 


“0. 004 etc. 
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'Compute Coefficients 


Compute 

da/dv 

Compute 

dM 4 /dv 

Compute 

dio/dv 

Compute 

de/dv 

Compute 

di/dv 

Compute 

dft/dv 
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fig. 4-3. FLOW chart of symbolic integrator 
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APPENDIX A J 

Effects of Earth's Oblateness 

A satellite, in an orbit of high eccentricity, spends 
a considerable portion of its orbital revolution quite far 
away from the earth (on the apogee side of the orbit). Never- 
theless, the effects of the oblateness of the earth are signi- 
ficant. The secular effects of the principal term, J 2 0» was 

the sole term considered. Note that J 30 a^d subsequent terms 

“3 ' 

are at least of order 10 compared to J 2 o 

It is well known that there are no secular variations 
due to J 2 o > over ^ unperturbed orbital period of the satellite, 
in tne semi-major axis, eccentricity and inclination. The 
variations in the longitude of nodes and the argument of perigee 
are found in many books on celestial mechanics. The variation 
for the last element M*, however, is not given elsewhere in the 
fora presented here, to the author's knowledge. The integra- 
tion, of this element is dealt with in detail. 



A1.2 


Ob la t oriels: Potential and Force 


For a spheroid (symmetry with respect to polar axis), 
the potential is given as 



R 


m 


J (— ) P (sin <5)) 
m=2 m r m 


(A.i.l) 


where 

r ' - distance from the center of mass of the earth 

to the satellite 

J - numerical coefficients 

m 

R = radius of the earth 

e 

P^(sin <5) = m-th Legendre Polynomial in the argument (sin 6) 
5 = latitude of the satellite 

and y = gravitational constant of the earth* 


The values of the first few constants J 

m 

are: 

—6 

J 2 - 1082.86 ± 0.1 x 10 

J 3 = “2. 45 ± 0.07 x 10' 6 

-6 

J 4 = -1.03 ± 0.2 x 10 


Consider the m-th term 

U « 
m 


K 

m 

r m+l 


P (sin S) 
m 


(A. 1.2) 


K 

m 


pJR 

m e 


where 



Then, the force due to this term is 


F = VU 


m 


m 

-K P 7(-“T) 
m m r IT,+ 1 


K + 
““ V(P ) 
m -H m 


(A. 1.3) 


= - 


, ( m + U . T 


_m+2 


(A. 1 . 4) 


^^ere 1 is the unit vector directed to the satellite* 


■ E *i 7T p m 


where x = unit vectors of the earth reference system (referred 

i 

to, elsev/here, as x a> y^, z q ) 


(P„) = 2x, P^ ■— (sin 6) 


m 


i m 3x. 

l 


(A.l .5) 


where 
Note that 


P* (z) = P (z) 
m dz m 


sin 5 = (1 x ) 
r 3 


Carrying out the algebra, Eq. (A, 8.5) reduces to 

' sin 6 * ' *3 


P = -P 


1 + P 


m m r r m r 


(A.l.. 6) 


Now, 


x = (x 1)1 + (x 1)1 + ( x 1)1- (A. 1.7) 

■a ^rr 3 n n 


where <l r . l t . l n ) are the instantaneous orbital axes and 
the parentheses indicate dot products. 



M.4 


Also, from orbital geometry, 


(x l r ) = 

3 r 

V - 

& I n ) - 

3 n 


sin 6 - sin i sin (co+v) 
sin i cos(u+v) 
cos i 


(A. 1.8) 


Substituting Eqs. (A. 1.7) and (A. 1.8) in Eq. (A. 8. 6) and 
simplifying , 




m 


p 

- —[sin i cos(u>+v)l_ + cos i 1 ] (A.i.9) 

r t n 


Substituting Eqs, (A, 1.4) and (A, 1.S) in Eq . (A. 1.3), 


E 


m 


= VU 


m 


K 

— “*[ (m-M)P .1 
r ra+2 m r 


P (sin i cos(w+v)l. + cos i 1 )] (A, 

m t n 


where the argument of the LP, and P , is sin 6. The prime 
denotes differentiation of P^ with respect to the argument. 
The components of along the instantaneous orbital axes can 
be determined by taking the appropriate dot products. Thus, 


F = 


(F 1 ) 

m r 


K <m+l) 

m ' ' , m +2 


m+2 


m 


m .m+2 . , , . 

F = (F 1 ) = - — A P sin i cos (u+v) 

? v m t m+2 m 


F = 


(F 1 ) 

m n 


m+2 * 

A P cos i 


m +2 


m 


(A, 1.11) 


J.10) 



where 


A = (1 + e cos v) 
and r ** p/A 

Secular Var iation in Q 

(A. 1.12) 

integrating with 

(A. 1.13) 
(A. 1.14) 

Secular Variation in qj 

^ ~ [-F + F - sin v] - (cos i) “ (A.l .15) 

dv ye 1 x 2 P dv 

Note that 

F - F cos v - F sin v 
x i 2 

Substituting from Eq. (A. 1.11) for n=2, and integrating with 
respect to v between the limits 0 and 2ir , 


dQ 


F sin(to+v) 


dv yae(sin i) 3 

* 

Substituting from Eq. (A. .11) for m=2, and 
respect to v between the limits 0 and 2r , . 

3K„ 


(An) 


ob 


~2 TT COS i 


VP 


Letting 


- K o R e 2 

K. - ~~2 ~ J 

PP 2 2 P 


(Afi) , = -3K ii (cos i) 
ob 



A1.6 


K,ir 

(Aoj) = - ~rr (18 siu 2 i - 12) - (cos i) (Aft) 

ob -4yp z 

Substituting for Afl from Eq. (A. 1.13) and for ^ from Eq. (A. 1.14), 
and simplifying 

(A o>) , = 6 Ktt(1 - 7 sin 2 i) (A. 1.16) 

ob .4 

Secular Variation in 


dM * = -2r 3 + 3 nt da 

dv 1/2 i 2 s dv 

yae 


1/2 

e 


(cos i 


dQ . dm 
dv' + d^> 


where n - mean motion of satellite 


and 


2 

da _ 2r a 
dv * ye 


( eF y + V 


(A. 1.17) 


This integration will be looked at term-by-term* 
Substituting from Eq, (A.iLll) for n=2. 


Term (1) “ - 


2r 3 

1/2 

yae 


F 

1 


2r ' 


yae 


1/2 


SKj, ^ 

P 


P 

2 


» -6Kc 1/2 P A 
2 

1/2 

= -3Ke A(3sin 2 i sin 2 (to+v) - 1) 


for P = P (sin 6) and sin 6= sin i sin(u+v) 
2 2 



Thus, on integration, 


Term (1) 





1 


- - sin 2 i- cos 2(w+v))dv 

2 


V 2 /3 

= - 6Ke it ( 2 


sin 2 i - 1) 


Next , 


Term (2) 


3 nt da 
2 a dv 


3 nr da 3 n(At) da 

“2 2a dv '2 a dv 


where 


and 


t = time measured from perigee 

+ c^) 

x “ period of satellite 
(At) — time measured from apogee 


Since the secular variation in a is zero, the 
ten. on the r I ,,ht-han,l side of lit,. (A.i.'l*) drops out 
** rat Ion. 

2 

Term (2) = 3n(At) ^ (cF y + F^) 

m IBeIIALL (eF S in v + F A) 
ye A 2 1 2 

F 

y 


(AJ .18) 


(A. 1.19) 


first 
on in Le- 


as 


= F sin v + F cos v 
1 ^ 


(A. .1.20) 



9nIC p - . „ 2 

. Term (2) =.“ (At) [-'(3 sin 2 i sin 2 (urt-v)-l) A sin v 


sin 2 i sin 2(m+v)] (A. 1.21) 


Integration of Eq. (A. 1.21) will be written out in detail, 
The generic form would be 


where 


/ (At)f (v)dv - (At) / f(v)dv - - /% (/ f(v)dv) 

C 1 A 


dv 


3/2 

P 1 11 

(At) « — ~ T = " “2 


A Z c/ 


(A. 1.22) 


Thus, (the limits, 0 and 2 tt , are not marked for convenience) 


- | / (At) A^sin v dv 


_ e r _,(A.t)A2. 
2 l ^ 3e ; 


3eC 


/ A dv] 


x ,, . x 2 tt 

6 (1 + <i) "15 


1 


(A. 1.23) 


3e 


sin 2 i sin 2 w /(At) A 2 sin v cos 2 v dv 


= —■ sin 2 i —■y— • / (At) A 2 (A 2 - 2A + l)sin v dv 


3e 


sm 


2 i sin 2 w ^- + - e o (-1 + 3e - 6e 2 ) 


30e' 


+ — V (2 + 3e 2 ) ] 
30e Cj 


(A. 1.24) 
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~ sin 2 i cos Z w /(At) A sin 3 v dv 


i~ sin 2 i cos 2 w / (At) A 2 sin v (l-cos 2 v)dv 


= sin 2 i cos 2 ^[- T - ^ -"^ (l-3e-4e 2 ) + — ~(-2+17e 2 )) (A. 1.25) 


30e 


30 b C ^ 


i— sin 2 i sin 2d) /( At)A 2 sin 2 v cos v dv 


= 0 as (At) is odd about v=ir 


sin 2 i 


sin 2cu /(At) A 3 (cos 2 v - sin 2 v)dv 


- sin 2 ! cos 2o) /(At) A 3 sin v cos v dv 


= - sin 2 i cos (l-4c) + <?e 2 -2) ] (A. 1.26) 

20 e 20 e Ci 


Putting Eqs . (A.. .23) to (A. .26) together, and simplifying, 

Term (2) =-^[777 (3 sin 2 i - 2)+ 7 (l+e) (1-3 sin 2 i sin 2 w)] (A. 1,27) 
N ' e 6 C 1 o 


Finally, 

Term (3) = -e ' (cos i (Afi) + (Aw)) 


1/2 


~ KTre ly,2 (sin 2 i - \) 


(A. 1.28) 
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Then 

(AM^) ^ = secular variation in due to J 2 o 
= Terra (1) + Term (2) + Term (3) 

Summing Eqs. (A. 1.18), (A. 1.27) and (A. 1.28) 

(AM^)^ = -3Kt 1/,2 rr(3 sin 2 i - 2) 

+ |ktt (3 sin 2 i - 2) 

+ |kttg 1 ^ 2 (3 sin 2 i - 2) 

+ | nT^ (1+e) (1-3 sin 2 i sin 2 m) 

Substituting -yr = e*^ 2 and nt * 2 tt , 

E 1 

(AM.) , = (1 + e) (1-3 sin 2 i sin 2 m) (A. 1.29) 

* ob e 

Features of AM,., 

There are several interesting features in the final ex- 
pression for the secular variation in M* over an orbital 

period. of the satellite, as given by Eq. (A. 8. 29). Firstly, 

2 

R J,R 

p e a e 

Since a is relatively invariant in high eccentricity orbits, 
Eq. (A. 1.30) shows that e (being of magnitude 0.1 ^ 0.2) acts 
as a small divisor in (AM*) b (as also in (A0) ob (Aco)^). 



Al.ll 


The effect is rather pronounced in (AM*) ^ because is pre- 
sent while, in (AQ) ob and (Au>) ob , c 2 is present. Thus small 
inaccuracies in e, the eccentricity, are magnified. 

Secondly, because of the presence of in the denomina- 
tor, the magnitude of (AM ft ) b is rather large. This is con- 
firmed by experimental results as shown in Chapter 2. 

Thirdly, the final expression for AM* arises exclusively 
from Terra (2) in the integration process (Eq. (A. 1.27)). The 
contributions of Term (1) (Eq, (A. 1.18)) and Term (3) 

(Eq. (A. 1.28)) cancel with. part of Eq. (A. 1.27). Thus, the 
"short-period” variations in the semi-major axis contribute 
to the secular variation of M*. 

Fourthly, the secular variation in M* is zero only when 


1-3 sin 2 i sin 1 '!) = 0 


Secular Variation of Roth's Sixth Element 
From Eq. (2.21) of Chapter 2, 

4 

—• = r , — [-Ft cos v + (1 + -)F 0 sin v] (A. 1.31) 

dv pe/pp 1 p *■ 

Substituting from Eq. (A. - .11) and after some manipulation 


dT 

dv 


Ke 


3/2 


ne 


i i 

[-3P. cos v -P sin i sin v cos(to+v)(l + -)] (A,. 1.32) 


Instead of transforming to E, the eccentric anomaly, as the in- 
dependent variable, this expression will be integrated with 
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respect to the true anomaly. The result in either case will 
be the same so far as the secular variation over an orbit is 
concerned. Detailing the integration term-by-term, the limits 
being assumed, 

3 

Term (1) = - /(| sin 2 i sin 2 (to+v) cos v - ^ cos v]dv (A. 1.33) 
= 0 (A. i. 33) 

Term (2) = - / - sin 2 i sin v sin 2(aH-v)dv 

= 0 (A. 1.34) 

Term (3) = - - sin 2 i / ^-[sin v sin 2(«+v)] 

Splitting Term (3) further, 

q T_rv \) 9 

Term (3.1) = sin 2m / — “ — cos'v dv = 0 

ei n ^\) 

Term (3.2) = -sin 2u / — “ — dv = 0 

„ sin 2 v , 

Term (3.3) = 2 cos 2w f — “ — cos v dv 

= 2 cos 2w[/ dv - / dv] 

= 2 cos 2ui [- I - ( ' A ~ 1 ' ) - dv 
e A 

- / (A 3 -3A 2 +3A-1)~J 

e 

Substituting 

, dv _ 2tt 
! t 


between the limits 
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and carrying out some simplifications, 

cos 2 m 

Term (3.3) = - 2ir 3 — (1 - t ) 

e 

Summing Terms (3.1), (3.2) and (3.3), 

Term (3) = ^ sin 2 i cos 2m (1-e ^ ) (A. 1.35) 

e 

Summing Eq. (A. 1.33) to Eq. (A.i.35) and substituting in the 
integration of Eq. (A. 8. 32), 

3/2 1/2 2 

(AT) t = 3 ttK ~ — r (X-t ' ) sin 2 i cos 2m (A. 1.36) 

ob ne 4 

This is the final expression for the secular variation of the 
time at epoch, due to J 20 J 2 s obtained from Roth s sixth element. 
Notice that e does not appear as a small divisor unlike in 
Eq. (A. 1.31). Thus, the magnitude of (AT) ob would not be large^. 



FIG. Aj. 


COMPARISON OF ROTH’S SIXTH “ELEMENT 
WITH OUR SIXTH ELEMENT. 

{OBLATENESS ONLY) 

(note: quantity plotted on the 

Y-AXIS IS THE MAGNITUDE 
OF THE ERROR.) 
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CHAPTER 5 

Singularity-free Methods, Using Regularization,, 
for Circular and Elliptic Orbits 

5.1 Introduction and Motivation 

Under this grant, singularity-free methods of orbit calculation, 
using regularization, and applicable to both the circular (e— 0) and 
elliptical orbits, have been developed and comparatively studied, 

in a Ph.D. thesis by S.K. Bhate^ 5 as the Prin- 

cipal Investigator's (M.L. Renard) advisee. For a much more detailed 
treatment, the reader should refer to Ref. [5-1]. Although the ori- 
ginal grant, in April 1968, had the title "Launch Window Analysis of 
Highly Eccentric Orbits", it became readily apparent that for missions 
such as that of IMP-H, methods of orbital and mission analysis for 
large circular orbits were required, which should be insensitive to 
the following singularities introduced by the choice of the standard 
"osculating elements" 

a, e, i, oo, n, T 0 (5.1-1) 

a) e = 0. The orbit is exactly circular. Since perigee 
strictly does not exist, the "argument of perigee" or 
"time of passage at perigee" lose their meaning. Mathe- 
matically, in the equations expressing the time-rate of 
change of the osculating elements, (5.1—1) > small or zero 
divisors "e" appear . 

b) i = 0° or 180°. The orbit is exactly equatorial, posi- 
grade or retrograde. Since the line of nodes strictly 
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does not exist, the "nodal line” is left undefined. 
Mathematically, small or zero divisors "sin i" appear. 

If the small, perturbing forces do not derive from a potential, 
or if such forces might have to be encompassed by the theory, it 
appears normal to write the equations expressing the time-variations 
of the elements in Gaussian form * The problem of developing a varia- 
tion of parameters scheme in which derivatives are expressed in terms 

r ^^ 2 1 f 5— 3 1 

of perturbing forces was attempted by A.M. Garafalo 1 , R.R. Newton 1 , 

[ 5—4 1 

CoJ. Cohen and E.C. Hubbard . In all of these studies, however, 

absence of any perturbations, which defeats the very purpose of the 

r 5-si 

method of variation of parameters. S. Pines presented the first 

"authentic" variation of parameters scheme, which used as osculating 
elements the position and velocity vector at some instant, time being 

and as the independent variable. Basically, the same method was used 

r, fil [5-7 ] [5-8, 5-9] 

later by P. Wong L , S, Herrick , E.', Pitkin *. Program NICE-T 

developed under this grant at C-MU , and described in 

Chapter 6 of this work, has been written based on Pitkin’s version 

of the variation of parameters, and will be compared to other methods 

presented in this chapter. 

In the following, non-singular elements such as the radius and 
velocity vector at some epoch are combined with the use of differen- 
tial transformation of the independent to result in the extremely 
simple form for the unperturbed equations of motion; those of the 
harmonic oscillator. Based on this unperturbed solution, a singular 
free method of variation of parameters can be developed. This is the 



5-3 


object of Section 5.2. In Section 5.3, a modification of Brouwer’s 
method of "perturbations in rectangular coordinates"^ 5 10 ^ which is 
applicable to circular orbits is developed, again starting from the 
unperturbed solution mentioned before. In Section 5.4, the perturb 
ing forces due to a third gravitational body (such as the sun or the 
moon) are expressed in terms of mixed Fourier-Chebychev series, for 
which series computational algorithms are derived which maintain a 
good accuracy by avoiding the problem of taking the differences of 
large, close numbers. This leads to the development of a theory 
which is semi— analytic , namely closed form integration is performed 
on series of the type indicated, the coefficients of which are 
numbers "valid" over one orbit of the perturbed body (the satellite). 

In Section 5.5, a numerical comparison is made between the results 
of the integration, in their various forms, of the system of differ- 
ential equations . Two "benchmark examples are considered, an orbit 
of large eccentricity (e = 0.936227) with a period of 4.45 days, and 
an essentially circular orbit (e = 0.8018212 x 10 ) having a period 

of 12.05 days (orbital radius = 35 mean Earth radii). 

For such "large circular" or "high eccentricity" orbits, it should 
be stressed that: 

— as compared to close— to— earth orbits, the oblateness effect 
and atmospheric drag, as perturbing forces, are most of the time, or 
even always, much smaller than those due to the gravitational pertur 
bations of the sun and the moon. 

— as compared to classical problems of Celestial Mechanics 
dealing with natural satellites: the perturbing forces considered in 
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these problems are small either because of the large distances involved, 
compared to the orbital semi-major axis (in the case of the moon) or 
on account of the small ratio of the masses (the largest for Jupiter 
but still < 10 -3 ). Typically, the ratio of the magnitude of the per- 
turbing force to that of the central force is largest in the case of 
Mars perturbed by Jupiter, and is .763 x 10 at most. 

To fix the ideas, the relative order of magnitudes of the per- 
turbing forces, with the central force as a norm, are, along large 
circular orbits: 

-27 

Atmospheric pressure < 10 

— 1 6 

Radiation pressure < 10 


Oblateness 


< 10 


.“5 


— ^ 

Perturbation due to the sun ^ 10 


Perturbation due to the moon ^ 10 


-2 


Now, in studying orbits, the main motivation might be: 


1) High accuracy computation of ephemerides : for this, numerical 

integration is well suited and can be carried out to a very high 
degree of precision. 

2) Determination of the evolution of the orbital elements, for mission 

analysis purposes: Here, since the requirement on accuracy is re- 

laxed, it might be allowed to linearly superpose perturbations, at 
least over one orbit of the satellite, in spite of the significant 

magnitude of some of the forces listed above. One might then 

[5-10 to 5-15] 


consider to use close-earth satellite theories 


Lunar 



5-5 


and planetary theories would have to be excluded if the perturbing 
function is expanded in powers of the inclination as a small 
parameter, and if the theory is therefore unsuitable for the large 
inclination commonly encountered for artificial satellites. A 
notable exception is Tisserand T s^ 5 theory for the computation 

of Pallas perturbed by Jupiter, in which the perturbing function 
is expanded in powers of both sin" and cos^ 


3) Orbit classification: Here the emphasis is on qualitatively 

classifying orbits, such as being able to say if they are of 
circulatory (line of apsides rotates monotonously) or oscillatory 
(line of apsides oscillates between limits) nature ^ ^ To this 

effect, some sort of development in series (Legendre Polynomials 
and Taylor series^ 5_17 ^ , Fourier series in M and M 18,5 19 ^ ) 
is considered and some vt main" contribution is analyzed to define a 
qualitative behavior of classes of orbits. 


The present chapter strive for the development of methods of rela- 
tively moderate accuracy and concentrates on this objective along the 
lines described in item 2). 

. As a last remark in this introduction, we should mention that we 

initially proposed, in 1969, that the study of large circular geocentric 

[ 5-201 

orbits be a part of the material studied under this grant , for 

orbits having a ratio of the orbital semi-major axis to the semi-major, 
axis of the moon up to about 7 ; , along the following lines of effort: 
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- Regularization (for instance, Kustaanheinp-S tiefel ! s trans- 
formation) 

- Closed form theories . 

These are indeed, the directing lines taken in this chapter* However, 
in view of the magnitude of the dynamical perturbations, the "analyti- 
cal" integration of the equations of motion was carried out with 
numerical coefficients (as opposed to literal coefficients) inside 
the computer program and its result evaluated to give the desired 
output. Hence the name "semi- analytic integration" was thought to be 
more appropriate. 


5*2 Unperturbed and perturbed two-body motion 

5*2.1 Development of the linear equations for the unperturbed problem 
In its classical form, the unperturbed two-body problem, referred 
to the center of the Earth, say, is described by the non-linear differ- 
ential equations 

r + - 0 (5.2-1) 

r 


in which r is the geocentric vector to the perturbed body, and y = k 

m + m ) or u ^ k 2 M , if m <<< M ... If a 

^ earth satellite'' ^ = earth satellite earth 

perturbing force F, assumed to be always of "small" magnitude compared 
to the magnitude of the central force, is present, (5.2-1) has r.h. side 


-j- 

F. 


Using the well-known Sundman’s transformation 

d .Id 
dx - ^ dt 


( 5 . 2 - 2 ) 
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leading to 


d_ = ^d_ . d£_ _ {v_ d_ /j± d_ 

dt r dx dt 2 r dx r dx 


y d 2 _ y dr d 
dx 2 r3 dx dx 


Equation (5.2-1) from which the consideration of collision orbits is 

. 2 - 1 

excluded (thus r ^ 0), will read after multiplication by r y , 


d 2 r 



- (—•) (%) + - = 0 
* dx dx r 


(5.2-3) 


which is still non-linear. Now, in two-body motion with an inverse- 
square attractive law of forces, and with c = r Ar = constant vector 
( f are derivatives with respect to time) 


r Ac -^(r Ac) 


-> -> t 

r-r 

v(x — r~ 


-> T 1 . 
- r — ) 

r 


->i -y 

/ r r . d / r \ 

= u( 7 - r - ) - v -£(-) 


Therefore, upon integration, 


eP - r A(r A r ) -p — 


(5.2-4) 


Let t. be the Laplace vector 


**>■ 

■f e ^ r 1 -*» x 

A = --P = -- -r A(r A r ) 

def V r y 


(5.2-5) 


It is a constant in the unperturbed motion. So is the energy inte- 
gral (a > 0 for elliptic orbits) 
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( 5 . 2 - 6 ) 


or equivalently the quantity 

a y y dt dt r a 


(5.2-7) 


Now A and a are rewritten in x variables as 


r 1 dr A *+ h dr. 
A r — A (r A —) 


r dx 


dx 


a = 


1_ dr_^ dr _ 2 
r 2 dx dx r 


From these expressions, otr + A can be expressed as 


*> -y -->■ ,■>* 

-+ ? 2r 1 ,dr dr dr dr -K r_ 

ar + A + — (— r • — r - ~r • — r)+ - 

r dx dx dx dx 


. 1 /dr -*\dr 

+ -y(— * r)-r 

r z ax dx 


_r + 1 dr dr 
' " r r dx dx 


Equation (5.2-3) is rewritten 


d 2 ? t 

— ry - ar = A 
dx 2 


(5.2-8) 


(5.2-9) 


( 5 . 2 - 10 ) 


(5.2-11) 
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Defining a new, auxiliary vector which is constant in the unperturbed 
motion. 


-+• ?■ 
B = aA 

def 


( 5 . 2 - 12 ) 


we obtain 


- a (r - 


(5.2-13) 


A variable z is defined as 


z = r - B = r -aA 
def 


(5.2-14) 


and the final system of uncoupled, constant coefficient system of 
linear differential equations reads 


d 2 z 1 
— T^T + “ z = 0 


(5.2-15) 


Its solution is, in terms of the z variable 


js . , r~ A z \ ■ ™ 

z = z 0 cos / H /a (— ) sm F~ 

/a dx o /a 


dz _ . 2£ , ydz N 2£ 

— = - / — sm r~ + (~T“) cos /— 
Hv v a Vet dx o v a 


(5 . 2-16) 


4 dz * , 

in which z„ , (—7-) are the initial conditions, given at x 
dx ° 


In terms 


■± .-*■ + ./— (— r) sin j— 

= B + (r 0 - B)cos y— +/ a dx ° * a 


dr r °~ B . . ,dr. £ 

— r sin 7- + (— ) cos 

dx / a /a dx o / a 


(5.2-17) 
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(5.2-17) is the solution of the unperturbed two-body problem in terms 

A 

of independent variable x. The time t corresponding to any x, is 
if t = t 0 for x = x 0 = 0, 


t 


x 


M d x 


J /y 


+ t c 


(5.2-18) 


Looking at Equation (5.2-17), it appears natural to introduce a new 

X 

variable /= . However, since a is not a constant in the perturbed 

v a 

motion, the introduction of the new variable "E" is done in differ- 
ential form: 

A 

dE = ¥ (5.2-19) 

ya 


In the unperturbed problem, considering that a is a constant of mo- 
tion, by integrating (5.2-19), we obtain 


E 



with E = 0 for x = 0 


( 5 . 2 - 20 ) 


whereas for the perturbed problem 
r E 

x = /a dE with x = 0 for E = 0 

0 

With this new variable. 


d dE d 

dx <*x dE 


d 2 I d_. 1 d_, 

d ~2 JZ dE^/a dE ; 

JL,d£ JL_ da d* 

" a dE 2 ~ 2a dE dE' 


( 5 . 2 - 21 ) 


(5.2-22) 
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In the unperturbed problem, “ is zero. Thus Equations (5.2-3), 

ah 

(5.2-11) and (5.2-15) become 


d 2 t _ldrdrr = 
dE 2 r dE dE r 


d 2 r 

dE 2 


+ r = B 


d 2 z 

dE 2 


+ z = 0 


(5.2-23) 


(5.2-24) 


(5.2-26) 


Therefore, in variable "E" and for the unperturbed problem. 


z = z 0 cos E + 


S) 

ME ' 0 


sin E 


dz -*■ „ , ,dz. „ 

— = - z D sin E + ( — ) cos E 

dE dE o 


(5.2-27) 


r = B + (r„ - B) cos E + (~) sin E 

ah O 


it = ^ Sin E + C0S E 


To obtain the differential equation satisfied by the scalar r, the 


identity 

2 ■+ -*• 

r = r-r 


is differentiated twice with respect to E 


r 


dr 

dE 


-*■ 

r • 


dr_ 

dE 


dnr 

dE 2 


-> 2 
,dr. 2 d r . 

(Ttr) = r 2 + ( 

dE dE 


dr 

dE 


2 
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From Equation (5.2-23) 



dr., dr 

dE^dE " ar 


d r 2. 
(— ) - 

W 


ar 


Substituting, and dividing by r, 


dE 2 


, l,dr s 2 

a + ; ( d5> 


From (5.2-9) 


a 


1_ dr 

0 

r dx 


dr 2 _ 1 dr dr 2 

dx ~ ^ = ar 2 dE * dE " ^ 


we obtain in Equation (5.2-29) 


(5.2-29) 


i ( d £) 2 

rW 


= ara + 2a = - r + 2a 


Hence 


dE 2 


+ r - a = 0 


or with z = r - a, 
def 


d 2 z 

dE 2 


+ z = 0 


the solutions of which are 


z = ^ sin E + c 2 cos E 


— = c 1 cos E - c 2 sin E 


(5.2-30) 


(5.2-31) 


(5.2-32) 
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and 


dr 

a + Cj sin E + c, cos E = a + (r 0 -a)cos E + (^) 0 s; *- n E 


dr _ _ 

— = ci cos E - C 2 sm E 


(5.2-33) 


with c^ c 2 constants equal to 


/dr N 

C 1 “ ^dE o 


C 2 


Note that the above equations (5.2-32), (5.2-33) will still hold true 

dx* 

for a perturbed problem, if a, r 0 and (r^) are osculating parameters 
of the orbit at any given time. 

The equation relating time to E, in the unperturbed problem, is 


t-t 0 = 


) 


/a dE with t = t e for E = 0 


or 


a -(. 2 - [E + (— - l)cos E + -(—•) (1 - cos E)] (5.2-34) 

t - t 0 = f— a adE« 

»V 


If the reference (t - t 0 ; E = 0) is the perigee of the unperturbed orbit. 


r 0 = a(l - eE) 
<§> - 0 


then. 
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Therefore, in (5.2-34), 

V 2 

a 

t - to = r [E - e sin E] 
vy 


(5.2-35) 


which is the classical form of Kepler's equation of time, E being the 
"eccentric anomaly". (5.2-34) this appears as a more general equa- 
tion, or generalized equation of time, with E playing the role of a 
generalized eccentric anomaly. 


5.2.2 Differential equations in the perturbed problem 

If a perturbing force F exists, Equation (5.2-3) written with x 
as the independent variable will read 


dfl . 1 dr d£ ? , I 2 

0 -V- A. \ ^ E 

dx Z r dx dx r y 


(5.2- 36) 


Similarly, with E as independent variable, and taking Equation (5.2-22) 


into account, 

d 2 r 1 dr .dr. r _ F 2 1 da dr 

dE 2 _ r dE ME' r y r 2a dE dE 


(5.2-37) 


or 


r 


dE" 


. * -> F 2 . 1 da dr 

+ r - B - — r a + — — — 
y 2a dE eE 


(5.2-38) 


5,2.3 Formulation of the variation of parameters method, with E as 
independent variable 
5. 2, 3.1 Variational equations 

Following Lagrange 1 s method of variation of constants, the solution 
to the perturbed oroblen is written as 
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r = B + (r 0 - B)cos E + (^r)^sin E 


(5.2- 39) 


4§ = -<r„ “ 5)sin E + (g) cos E 


dr. 


dE 


dE 


If at the instant considered, i.e. for the value of E considered, 
the perturbations were removed, vectors B, r Q and (~j^r) would be 
constants of the motion in the unperturbed motion that would ensue. 
However, due to the continuing action of F, B, r 0 and (^g) will be func- 
tions of E, the variations of which are now determined. 


Taking the derivative of r in Equation (5.2-39), 


4i = #U - cos E) + cos E + )sin E"(? c -B)sin E 


dr 


dE dE 


dE 


dE dE 


dr 

+ (^) o cos E 


(5.2- 40) 


Subtracting the second equation of (5.2-35) from (5.2-40), we obtain 


{§<i - « + zr E + ^«S ) - >5ln E ‘ 0 


(5.2-41) 


d^r -4 ”> 

The l*h. side of Equation (5*2-36) is equal to ~o + r “ ^5 thus in the 

dE'" 

perturbed problem 


^i + ? 

dE 


s - % sln E - if 1 sin E E 


F o , 1 da dr 
= U r a + 2a dE dE 
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cT? A A* 

In terms of and rewritten 


dr o „ d /( ,dr. . _ F 2 djB . „ 

-sin E — + cos E ^ r a “ d E Sin E 


1 da dr 
2a dE dE 


( 


Solving (5.2-41) and (5.2-42) for the six-vector 


dE 


4- 

( dr } 
ME ; 0 


J 


(5.2-42) 


, yields 


d_ 

dE 


(— ) 
MB' 1 


r 


cos E 
sin E 


-sin E 
cos E 


x 


dB /n \ 

- ai a ' cos E) 

dB . ^ , ,E\ 2 ,1 da dr 

- dE S1 ” E + ( U )r 3 + 


<5. 


2a dE dE 


or explicitly 


d? J 2 „ 1 da dr . „ , dB „■> 

= "ii rasinE '^Ii sxnE di ( ] 


dE 


^r((7§)o> " 7* Z * cos E + h It ft cos E - dl S±n E 


dE dE 


2a dE dE 


(5.2-44) 


This system of differential equations of sixth order defines the varia- 
tion of parameters scheme, describing the variation with E of the oscu- 
lating elements r 0 , (”^) ° , integrals of the unperturbed motion. 


da dB 

To be complete, there remains to compute ”, — , 


da da dt_ 

d?" dt dE 


2-43) 
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Now, by definition 


1 _ 1 dr dr 2 

a y dt dt r 


1 da 2 dr _ d 2 r 2 dr 


. 2 dt 


V dt dt 2 r 2 dt 


2 dr r . 2 dr 

- — ' “ 7- > + 72 IT 


y dt 


. dr 2 dr 2 dr 

y dt r 2 dt r 2 dt 


da 

dt 


y dt 


Thus 


da 

dE 


2s!(|. i£) 

y v dE ; 


To obtain 3- , since by definition 
dt 


-> -> 

f r 1 dr , dr N 

A -7— A(r A 77) 

r y dt ; dt 


B = a A 


dB 

dt 


da ~T 

d? A 


+ a 


dA 

dt 


dA 

dt 


1 dr 


r 


dt 


dr 

dt 


1 d z r 


^ dt' 


a , dr. 1 dr . 

Hr A ^ ) 7dT A(r 


i; 

A--) 


d 2 - 


dt' 


( 5 . 2 - 45 ) 


( 5 . 2 - 46 ) 



5-18 


df£ 

dt 2 


is substituted from (5.2-1) and triple products are expanded, 


dA _ 1 djr _ r d_r r_ dr _ 1 dr 1 .-K dr 

dt r dt r 2 dt r 2 dt r dt y^ r ' 


'dt 


l/t dr,->- 1 ,dr +.->■ l.-»- dr.£ 

- (F • — )r-~-(— * F)r + - (r — )F 

y v dt y v dt y dt 


Finally, 


-> 

dA 

dt 


l r/ ± + \dr , ,-*■ dr.:± 

" ;< (F ' r) dT + (r - 3F )f 


)— 

dt 


2 (?- §)r] 


(5.2-47) 


_ dB 

To compute -j^r , 


dB 

dE 


,F dr, 
'y dt' 


2a(r ‘ ~)<aA 


, ,F ~>vdr , dr. F 

r) + a(- • r)-r— + a(r 3 — )- 
y dt dt y 


(5.2-48) 


and 


dB 

dE 


*<?. If) -?> + .<* 


(5.2-49) 


Equation (5.2-43) together with Equations (5.2-46) and (5.2-49) 

are a complete formulation of the method of variation of parameters in 
-+ dr 

elements r D , (— ) , with E as independent variable 
cm 0 


5. 2. 3. 2 Equation of time 


Using equation 

t - to 



a; dE 


(t 


t Q for E = 0) 


(5.2-50) 
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in which necessarily r is given by 

dr 

r = a + (r e - a)cos E + (— X sin E 
as in (5.2-33), but this time with the elements a, r 0 , (^r^ 
functions of E. Substituting in (5.2- 5Q , 


fE 


t - t 0 = 


^ [a + r 0 ~a)cos E + (^—) 0 sin E]dE 


— 3 /? 

v'y (t-t 0 ) = a E - 


E 4=r(a /,2 )dE +/a(r 0 -a)sin E 
dh 


f E 


sin E -^—{/a(r 0 -a) } dE - 


/a(||) o cos e) 


J 


1 cos EdE 


or 


t-to = 


3/2 

IE + - l)sin E + !^f) e (l - cos E)] + 


' E M 

/u 


dE 




(5.2-51) 


with 

f (E) = - E |^(a 3 ^ 1 2 ) - sin E ^[/a(r 0 -a)] + cos E~[/ a (||) 0 ] (5.2-52) 

V 

Written in the form of (5.2-51), and by comparison with Equation ( 5 . 2 - 34 ) 
for the unperturbed motion, it can be seen that the "change in time" is 
given by the last two terms 
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/E 

t " t — f (E) 

perturbed unperturbed 7 ==--*- dE + 

i v U 

0 

/Z<4fU). - /I«J§). (5.2- 5J 

77 

By observation of Equation (5.2-53), it can be seen that in order to 

know the perturbed time, one more integration in E is needed (the 

order of the system being seven). This feature is equivalent to that 

of a double integration, a step which cannot be avoided at some point 

[ 5 - 1 

as commented upon by J. Kovalevsky 


5. 2. 3. 3 An approximate variation of parameters scheme: a approximately 
constant 

If the nature of forces is such that along the perturbed orbit 


da 

dt 



-+ 



(5.2-54) 


is approximately zero (i.e. always much smaller than the time-rate 
of the other elements), or if, in the absence of a priori knowledge of 
such smallness of , numerical experiments have shown that such was 

the case in the problem at hand, a simplified scheme can be developed 
along the lines described in what follows. Extreme caution has to 
be exercised, however, in making sure that (5.2-54) holds sufficiently 
well. 

Let E be defined here as in the unperturbed problem: 


A 

X 
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Thus 


dx rr 1 

dE /a + 2 


= da 
/a dE 


= /a (1 + 


1 E da . 

2 a dE ; 


d_ 

A 

dx 




d_ 

dE 


The simplification made is that j 2E is always very small compared 

to unity, thus 


dE _ A 
dx 1/3 

With this simplification, i.-e. ^ neglected in the differential 
transformation, the equations for the perturbed motion as given in 
Section 5. 2, 3.1 simply reduces to 


if? 

at 2 


+ r 


-> F 

- B = - 

U 


r 2 a 


and the Equation (5.2- /+ 4) for the variation of the parameters 


become 



dE 


.dB 

(1 - cos E)^g 


(- r 2 a) sin E 

y 


(5.2-55) 


IS !#■>.! =" sinE # + l r2a cosE 

The relation between time and E is given by (5.2- 3^. 


5.2.4 Formulation of the method of variation of parameters, with 
time as independent variable. 

The solutions developed for the unperturbed or perturbed motion, 
with E as independent variable, are now transformed to the case where 
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t (time) is the independent variable. 


5.2.4 Unperturbed motion 
One wishes to obtain 
r(t) = f(t) r 0 + g(t)(^-) o 


dr 

dt 


= f r e + 


.dr. 


«£>• 


( 5 . 2 - 556 ) 


* designating now derivatives with respect to t. Since, from 
Equation (5.2-38), ^ 

r = B + (r 0 - B) cos E 4- sin E 


Now 


B - constant in unperturbed motion 
r ar ,dr. . a , dr.dr 

‘ a r - r <dt> + ; (r dt'dt 


+ .1 2 1. ar .dr. dr 

, ar( - — <— >— 


•> a . . ar dr. .dr. 

- r .a -— ) +-^ 

Using the relation dt = dE —=~ , 

vu 


r = r e {1 - - (1 - cos E) } + 

r o 


7=r {Sin E + % ( dt ^ . a-cos E)}(§)„ 


Therefore, in Equations (5.2-56), 

f = 1 - — (1 - cos E) 
r 0 

r / a /a H r 

g = -Jp (sin E + -jp (p 0 (l - cos E) 


(5.2-57 ) 
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Computing the time-derivatives: 


df a . 1 

__ — sin E - — f — 

dt r Q r /a 




f * I- 1 ’ • n 

- - — sin E 

r 


r 0 /a , „ , J a ,dr N . „./y 

g = — _ — (cos E + -== (— ) sin E)— 1 

/ u /y dt ° r /a 


= T 2 - (cos E + ~ (~) 0 sin E) 
r /y dt 


To summarize. 


* / y /a . 

f = - — sin E 

r r 0 


g = ^(cos E +^5 (~) 0 Sin E) 
° x /p dt 


(5.2-58) 


From these expressions, an algorithm can be implemented, which for a 
given value of t, will require, to obtain the corresponding value of 
E, to solve the transcendental generalized Ke>l er equation (5.2-34) by 
some numerical method. 


5. 2. 4. 2 Perturbed motion 


Again the perturbed solution is written in the form 


r(t) = f r c + g(^)o 


* „“>■ 
dr _ -> , * .dr. 

dt = f r ° + g( dt )o 


(5.2-59 ) 
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* » 


where r 0 , C^)o are now functions of time, and f, g, f, g are func- 
tions given in 5.4.2. 

y T* 

From Equations (5.2-44) and dt = dE > and Equation (5.2— 49 ) » 

, „ _ , _ . _ dB da 

(5.2-45 ) for ^ , JF 


dxj 

dt 


= - - r/a sin E - 


—y 

1 da dr » « . d B ✓ t->\ 

— — : — - r — sm E + —(1-cos E) 
2/ya dt dt dt 


(5.2-60) 


,F . +,dr 


f^^|(al4) + a(M)| + a(r^ 


da = . d£ 

dt u dt 


Finally, (^~) 0 is needed. From Equation (5.2-44), 


d_r /djf) i _ d_r /dr^ F -. - a . 

dE L ME ; ° J dE U dt ; ° /JT 


d r ,dr. , r r c a , /d r x r dr 0 / a , r n da n 

= “U— )«] ~TT^ + 4t ; ° L dE /y 2/a~ dt J 


dt dt 


r r ° a r^L/ /dr. , f dr, A dr^ + 1_ da. , 

u E dt U dt ; ° J dt 0 r 0 dt + 2a dt ;j 


and 






(5.2-61) 


/I dr o . 1 da . /dr . 
4a dt 2a dt M dt ;o 


Given the form (5.2-59) to the solution of the perturbed problem, Equa' 
tions (5.2-60) and (5.2-61) define the equations for the variations 
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of parameters r 0 , (“) 0 with time as the independent variable. 

This form of the method of variation of parameters is implemented 
in program NICE-T, described in Chapter 6. 

5.2.5 Brief Comments 

To conclude this section, a few comments are made on other 

similar approaches taken by other workers. In his "Theory of Orbits , 

Szebehely gives a treatment of regularization, with an extensive 

bibliography. The attention is focused on the restricted three-body 

problem, two-body problem and the collision orbits (one-dimensional 

problem). In this, the Levi-Civita transformation, or the use of 

complex variables, is possible. A recent extension of the Levi- 

Cxvita transformation to 4— dimensional space (its extension to three- 

dimensional space not being possible) , by Kustaaniheimo and Stiefel 

[ 5 - 22 ] 

(KS transformation) is used by Stiefeland Scheifele , with the 

same independent variable as was used here. A disadvantage of using 
the KS transformation is that the degree of the differential equa- 
tion increases to ten , whereas the present treatment uses only six 
differential equations, plus one additional one to go back to physi- 
cal time, in a derivation which is thought to be simple and straight- 
forward. The method of variation of parameters developed by Stiefel 

r p o i 

and Scheifele 1 J las ten parameters, with the equations of 
constraints used as numerical checks during integration. Furthermore, 
and this is of particular importance in view of our goal to develop 

methods suitable also for circular orbits, Stiefeland Scheifele's 

[ 5 - 22 ] 


parameters 


make reference to perigee and as such are not suitable 
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[ 5 - 3 ] 

for circular orbits* Finally 9 as compared to Pitkin s develop 

ment of the method of variation of parameters with time as the inde- 
pendent variable using the perturbative operator technique, the 
method developed here is thought to be much more straightforward. 

5.3 A Modification of Brouwer's Method of Perturbations in Rectangular 
Coordinates, Applicable to Circular Orbits. 

5.3.1 Introduction 

In 1944, D. Brouwer ^ 5 published a method of calculating per- 
turbations in rectangular coordinates, which is described in Brouwer 

„r5-2if] . . 

and Clemence "Methods of Celestial Mechanics as being apart 

from Hansen's method, the only other one that "need to be considered 

seriously for application where the numerical values of the elements 

are used from the start, and where a precision compatible to that of 

. [ 5 - 25 ] 

observation is desired." The method was applied by M.S. Davis 
to compute the motion of Eunicke (first order). Recently, S.A. 
Hanid^ 5-26 ^ developed a second-order planetary theory using this method. 

In Brouwer's method, two main parts exist: the first one is to 
set up the differential equations of motion, the second one to integrate 
a suitable representation of the perturbing potential or forces, in 
expressions which are easily integrable. 

In order to assess Brouwer's method as regards to formulation of 
the differential equation, a comparison is made in Chapter 6 between 
Brouwer's approach and the classical variation of parameters method 
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method using the same package for computing the forces and identica l, 
algorithms for the integration. On this very limited sample, it appears 
that a slight advantage might exist for Brouwer's method. 


5.3.2 Differential equations in a form valid also for circular orbits 
Since Brouwer's method was using a reference orbit given in 
Delaunay's variables, it cannot be used for circular orbits. A new 
method is developed on the basis of the equations of Section 5.2, 
namely the use of elements r 0 , r 0 with E as the independent variable. 
In Section 5.2, the variable 
z = r - B 

was introduced. In the unperturbed case, it satisfied the equation 


dE 


with subscript reminding us that it is the solution in the case 

where perturbations are removed. With the same notation, we can de 


fine 


So 


r 0 


~y 

z 


= "B" on the two-body reference orbit (a constant) 

def 


= r on two-body reference orbit, at E 
def 


= r - B 0 , r taken along the perturbed orbit 
def 


This gives in Equation (5.2-38 )» from which (5.3-1) is subtracted, 

da yF dr \ n 9 j ^ 1. n 4 J 


and in which = (“ * ^r) 2a 2 is substituted, 

- z° l + (z - t 0 ) = B - B 0 + | r 2 a + (- 
dE y P 


dr, dr 
dE )a dE 
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Let 



the vector difference between 
actual and reference orbit 


<5 B =B - B, 


Thus 


2 -* 


2 (<5z) + 6z = fiS + - r 2 a + (- • tt) 
dE 2 V V flK 


F dr. dr 
dE J a dE 


(5.3-3) 


The solution to (5.3-1) can be written as in (5.2-27), with 5 
introduced to avoid confusion in the notations. 


->- -y 

z 0 - C sin E + D cos E 


,dz + -*■ 

(■j^)o = C cos E - D sin E 


(5.3-4) 


The solution to Equation (5.3-3) in its homogeneous form (r.h. side = 0) 
is of the same form as (5.3-4), thus 


dz = J K. Ifr 2 
i=l l 3C, 


(5.3-5) 


in which the K.'s are constant coefficients, and (i = 1,2,3) are 

the projections of ^ on axes X, Y, Z; C^(i = 4,5,6) are the projections 
—y 

of D on the same axes. 

In the perturbed case, one requires that (5.3-5) still be the solu- 
tion to (5.3-3), but now with the 1<_ ' s being 6 unknown functions of E, 


which can be determined with the additional requirement that the actual 
and osculating velocity be the same. 
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Thus, 


d /x'N 
dE (SZ) 


E K — + E 2-^2. 

" i=l :L dE o C . i=l dE 8 C 4 


(5.3-6) 


term appearing 
in unperturbed pro- 
blem 


and it is required that 

l = o 

i=l dE 3C ± 

Further differentiation of (5.3-6) gives 


6 dK-s a 
£ / 1\ ° 


,dz< 


,d 2 t 


5* * s > ■ j, <dT> fc; 'as > + A K i ir t 


Replacing 0 
dE^ 


, 2-> 

d z 0 , ^ ~t 
by -z 0 , 


_ l , dK i>3 ,dz^ , y K 8 /_ t n 

' dE 2^ z) i=l dE 9 C ± ^dE * i=l K i SC. ( o) 


or using (5.3-5) 


il 

dE 2 


d^ _ 6 dK^ a /dz Q . 

?2 (6z) + 5z i £ 1 dE aChME ' 


Thus we have six equations for the 


dKi 

dE - 


£ dK i 3 / dz n . _ 

i=l dE 3C. ME ; 


' G x G y G Z >T 


(5.3-8) 


in which 

-> 

G 


F 2 , ,F dr. dr £ 

, = - r a + (- * tf) a — + o B 

def y y dE dE 


(5.3-9) 
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From Equation (5*3-4), if C has components (C^ C 2 , C 3 ) and D 
ponents (C^, C^, C^), the following holds 


3C. 

1 


<z 0 ,j) = sin E (i, j=l, 2, 3) 


= <5 cos E (i = 4, 5, 6; j = 1, 2, 3) 

(i-3) 


j? ( — z °- ') 

ac. me i 

i J 


Mz o y 

ac. ME •' j 
1 J 


6 J cos E (i, j = 1, 2 S 3) 
i 


-5. sin E (i = 4, 5, 6; j = 1, 2, 3) 
i~3 


In matrix form 


sin E 
0 
0 

cos E 
0 
0 


0 

sin E 
0 
0 

cos E 
0 


0 

0 

sin E 
0 
0 

cos E 


cos E 
0 
0 

-sin E 
0 
0 


0 

cos E 
0 
0 

-sin E 
0 


0 
0 

cos E 
0 
0 

-sin E 


\ 

0 


0 


0 


G x 


G y 

J 

G z 

J 


com- 


CS. 3-9) 



After inversion 


'V 


sin E 

0 

0 

cos E 

0 

0 


0 

*2 


0 

sin E 

0 

0 

cos E 

0 


0 

K 3 

— 

0 

0 

sin E 

0 

0 

cos E 


0 

*4 


cos E 

0 

0 

-sin E 

0 

0 


G 

X 

K 5 


0 

cos E 

0 

0 

-sin E 

0 


G y 

* 

CD 


0 

s 

0 

cos E 

0 

0 

-sin E 

J 


G 

zJ 


dKi 

dE~ 

dK 2 

dE 

dK 3 

dE - 

dK 4 

dE 

dK 5 

dE 


= cos E G 

x 

~ cos E G 

y 

- cos E G 

z 

= -sin E G x 

- -sin E G 

y 


(5.3-10) 


dK 6 

dE 


= -sin E G 


To go back to physical time, t, use can be made of Equation (5.2-48)* 
Alternatively , using 
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t - t 0 = 


_1 

vTT 


(ry^a - ro^o)^^ + y=r 


J £ 


dE 


(5.3-11) 


r 0 having here the meaning of r, along the reference orbit s at E. But 
it is known that, if ( ) 0 designates the value of the quantity between 
parenthesis at E = 0 , 

da 

r 0 = a c + ((r 0 ) o - a c ) cos E + ((— )o)o sin E 
Let, in (5.3-11) 

dK 7 _ 

— — ,= r/a - r„/a 0 (5.3-12) 

dE def 


Then the equation of time becomes 

1 


t - to = ~ 7 = 


G J 


dK? dE + {a 0 E + ( (r 0 ) o “ a e ) sin E 


dE 


y 


(5.3-13) 


A A 


+ ((—)o)o (1-COS E)> 


[ 5- 1 ] 

An algorithm can be developed on the basis of the above formulae 
and is implemented in program BROUWER-E described in Chapter 6. 


5.4 Semi-Analytic Integration Method: Mixed Fourier-Chebyshev Series 
5.4.1 Introduction 

Given a system of differential equations describing the rates of change 
of the parameters., the first step in Picard's iteration scheme will 
consist in substituting in the r.h. side of the equation the solution 
to the unperturbed problem, and proceeding to analytically integrate 
this r.h. side to obtain the changes in the parameters over a suit- 
ably selected interval. 

For the large circular orbits considered here, it has been found 
that results of sufficient accuracy (in a sense to be precised later) 
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are obtained over a range (0, 2ir) for E. For larger intervals of E, one 
proceeds to the integration over (0, 2 tt) , "updates” the elements by add- 
ing their changes over the interval to the initial values, and so on. 

The r.h. sides in the differential systems, in order to be able to 
integrate these numerically, are to be expressed in a series of suit- 
able analytic functions, the coefficients of which are numerical , intro 
duced from the start using the initial conditions. Such analytic 
functions are oftentimes double Fourier Series in the mean anomalies 
of perturbed and perturbing bodies. However, in the case of artificial 
satellites undergoing strong perturbations, it is no longer true that 
the elements will not "change" too much over a period of the perturb- 
ing body. From that viewpoint, the periodicity in M* , say (mean anomaly 
of the perturbing body), for constant a, e,..., has been destroyed, and 
it appears perfectly reasonable to use non-periodic approximation func- 
tions, valid over a suitable interval in E, to represent the motion 

of the perturbing body. Chebyshev's polynomials have been used in 

[5-27] 

planetary theory by Carpenter 

The analytic series chosen here to represent the terms in the r.h. 
side are mixed Fourier-Chebyshev series: the Fourier part accounts 
for the motion of the satellite, and the Chebyshev part, having an 
argument which has a linear relation with the variable in the Fourier 
Series, represents the motion of the perturbing body. 

5.4.2 Development of mixed Fourier-Chebyshev Series for the derivatives 

The point of departure is the system of equations for the variation 
d* 

of parameters r Q 5 with E as independent variable, written in 
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(5.2-44), which are repeated here 


dr c F o . 1 da dr . _ , , n 

= - - r^a sm E - — “ ” sxn E + (1 

y 2a dE dE 


„.dB 
c° S E)^ 


-> -4 ~4 ~4 

d //dr^ N F 2 r . 1 da dr dB 

dE dE 0 " y r a COS E + 2a dE dE COS E Sln E dE 


(5.4-1) 


4- -V’ 4- -4- ~4- 4* 

dB » ,F dr. , f 'N . ,F -K dr F, dr. 

dE ■ 2a ( ; • dE> (aA - r) + % • r)a dE + J< r di )a 


da r\ *• 

dE * 2a V 


n 4- 4- 

2 F dr 


dE 


with A (in this first order scheme) being a constant vector depending 
on the initial conditions. 

In order to determine what calculations are involved here in the 
development in series of the r.h. sides of (5.4-1), we shall, as 
announced in the beginning of this chapter, consider the analysis 
is limited to gravitational perturbing forces due to other bodies 


(moon and sun) : 


-4* 

F 


p=moon 5 

sun 


-4 "> 

U r t 

c/x-V 

M r 

IP 


(5.4-2) 


—4 -4 

Standard notations are used: r , r are the geocentric vectors to the 

perturbing body and satellite, respectively, and r is the magnitude 
of vector Ir -r I . 

p 

First of all, the r.h. sides of (5.4-1) involve scalar r and 

-4 

-4 dr 

vectors r, B, which areto he taken as. those of the undisturbed 

dE 


motion, namely 
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Substitution of (5.4-3) into (5.4-1) permits to express the r.h. 

F 

sides as Fourier Series in E, multiplied by - » In Equation (5,4-2) 
it can be seen that for given r p , - could be expressed as a Fourier 
series in E; r could be developed as a Fourier series in the mean 
anomaly of the perturbing body. However, here, these series lose their 
validity rapidly, if the range of E ( for which (5.4-3) is adopted with 
values of K§)o, B at E = 0, exceeds (0, 2ir). It was decided to 

adopt Chebychev series for representing r , by the method of special 
valpes^ 5-28 1 With a judicious choice of the argument of the Cheby- 
chev series, the integration of the resulting mixed series could also 
be simplified. 

Let x be an argument linearly related to E as folloitfs 

E — (x + 1) it (5.4-4) 

Thus, x has range (-1, +1) when E varies over (0, 2 tt). The ephemeris 
time, t, corresponding to E, is obtained from Equation (5.2-34) 

t - to = [E +(f JL " Decs E + i(~) 0 (l - cos E)] (5.4-5) 

r p corresponding to this t (or E) can be determined by any suitable 
means (epheme rides, tapes, tables) with a suitable interpolation routine 
(here: fourth order central differences). The numerical coefficients 
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of the expansion of r in Chebychev's polynomials follows the method 

[ 5 - 28 ] 

given by Fox and Parker 

Mixed Fourier-Chebychev Series will be expressions encountered 
in the r.h. side of (5.4-1) having the canonical form 


J K 

max max 

f - Z { E Cl (J,K)T .(x) }cos (J-l)E 

J “ 1 K— 1 

(5 . 4-6) 

dmax Kmax 

l { I S1(J,K)T (x)} sin JE 
J=1 K=1 K_1 


The following properties and calculations are derived in S.K. 
Bhate's thesis^" 5 1 ^ , to whom the reader should refer for a more de- 
tailed proofs: 


I) If f^, f^ are series such as in (5.4-6), so are 

f + f 
r l 2 

cf^ (c scalar 

gf^ (g Fourier Series in E) 


hf^ (h Chebychev series in x) 


£ 1 x £ 2 


(From this results that when the proper substitutions and operations 

are carried out, all r.h. sides in system (5.4-3) will be mixed Fourier- 

F 

Chebychev* s Series, provided the components of - themselves can be ex- 
panded in such series ♦ ) 


r 3 

-n J 


II) D 1 = H^-) and D = D -1 (r < r ) 

1 def r lp p def 1 p 
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are mixed Fourier-Chebychev series , 


Indeed, let p d | f — < 1, 6 = angle between r p and r 


r 3 
ip 


_i 


(r 2 + r p 2 - 2rr p cos 6) 2 (1 + p 2_ 2p 2 ±n _) 3 /2 


*0 o 

if a ,= e- 3 . Thus, after Taylor's series expansion about xi= p = 0, 

def 

, o-l A 

and x = p =0, 


/— = i — - — - t" — 

r 1 p (1 - pCF) 3 2 (1 - po X ) 3 2 


, ( f (an+l) l n n ( - gatlH ■> „-■>) 

n=0 2 2n n! n! m "° 2 2m „!»! 


“ “ , , n+m n-m 

E E A A p a 
n=0 m=0 n m 


with = ^k + ' ^ ' ! ' ^ k = n,m ^ 
^ 2 2k k!k! 


(5.4-7) 


Now, in Equation (5.4-7), the double summation is effected along lines, 
in an "n,m" plane, of constant r^ = n+m or r^ = n-m. (5.4-7) simplifies 
to 

r t)3 « 2 2n ?/"?/.. v urtri, ri 

(_E_) = j; Ap + E ( E (A A ) p ^ ) o 

v r J m=0 n r^l n=0 m m+r 




(m=n) 


+ S (2 (A A )p 2n+r 2)o -r 2 
r 2 =l n=o n n+r 2 



From the above expression for A^» it is easy to see that 


A = 1 


\+l (1 + k+l )A k 


and one can calculate, for given ra max > 
m 

a ° def m=0 m p 


m 


max 


a„ = 2 S A A . 0 p 2m integer > 0 

£ def m= 0 m m+£ 

r 3 

Thus, (—^—) can be rewritten 
r iP 


(_E_> 3 = a 0 + r | 1 a r p r cos (r0) 

r iP 

.[5-1] 


(5.4-8) 


(5.4-9) 


It can be shown 1 ^ J that the sum (5.4-9), truncated to N. terms, can 
be computed without running into the problem of subtracting two almost 
equal numbers of large magnitude, in the following manner* 

b N= 2 = ° 
b N' 1 = ° 

b = 2p cos 6 b -p 2 b „ + a (r = N, N-1,...,0) (5.4-10) 

r r+l r+z r 


SUK = bo ~ b x p cos 9 = | (bo - p 2 b 2 + a 0 ) 


As is evident from (5.4-8) and (5.4-10), the mixed Fourier-Chebychev 
series for p 2m and p cos 0 are needed. These are given in Appendix A- 2-1 
and A- 2 .2. Once these have been obtained, it has been proved that D 
and D as defined above, are mixed Fourier-Chebychev series. 
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TT 

III) The components of - can be expressed as mixed Fourier-Chebychev 

series. Indeed, with the above definitions, and assuming that “ has 

P 

been developed in Chebychev series by the method of special values as 
in Reference [5-28] , 


h 

y r. 


p p 


t , ii. [D.?_ - D,?l - ^ - B,r] 

y 


a 2 v r 


->• 

r 

P P 


p * ‘ P 

From what precedes, the component terms between brackets should be 

a 3 

Fourier-Chebychev series, and so is ( — ) . Thus, so are the components 

of - . 

U 

The other terms in the r.h. side of (5.4-1), these not involving 

-± 

- alone, should also be shown to the series of some type. Indeed, 


da 

dE 


„ ,F dr. _p 2 ,a , d -> uj. n 1 

= 2a ^ - (— ) [Djc n ' - V dE ] 


dr 


-*■ dr 


dE 


par 


P P 


dE 


is such a series. Indeed, r^ and r = are, from (5.4-3), the 

mixed Fourier-Chebychev series 


? P • # ' - - *» «*■ E + E 

r ' S" <E ' ( dE ) “ )cos E+ ZE 

-B‘C sin E + |[(||)«* (ff)o " C*C]sin 2E 


d| 

AZs0 5 dE 


41 - -?>s + (? • aa - ?) + e- 3>s> 


-> K dr 


dr, 




dE " p v r p / tL "p''p _/ dE ' '‘p dE 


dE p 


2 dr , dr 


-V r « + (r 


dE 


)(2B - r)} 


will be a series of same kind provided 
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r *r = (r -B) + (r • Qcos E + (r • (^) e )sin E 
P P P p dE 

is of the same kind, which holds true. 

In summary, the equations (5.4-1) for the variation of t 09 (—) 0 
can be written with each r.h. side in the form of a mixed Chebychev- 
Fourier series. The latter can be integrated analytically to yield 
the changes of the elements over an interval (0 5 2tt) for E. 

5.4.3 Integration of mixed Fourier-Chebychev series 

Ref. [5-1] develops in great detail the algorithms needed for 
the integration, with respect to E, of the mixed Fourier-Chebychev 
series present in the r.h. side of Equation (5.4-1). These algorithms 
are briefly listed in Appendix A- 23. 

The following should be noted 

- integration is to first order of Picard’s iteration. 

r 3 

- the series (— £— ) is suitably truncated 

r iP 

“ the motion of the perturbing bodies are represented by 
finite Chebychev series. 

- in series multiplication, truncation is effected without 
loss of accuracy 

This being said, no other approximation is introduced, since the inte- 
gration of the terms being kept is rigorous. 

A computer program, based on the above technique, is developed 
and has been tested for close to zero as well as for large eccentri- 
city orbits. It needs some more work, however, to incorporate time- 


saving shortcuts. 
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5.5 Conclusions 

In this chapter, methods of deriving variational equations for the 
elements of elliptic or circular orbits which are totally singularity- 
free (in the absence of collisions) have been developed. They should 
be particularly useful for the study of nearly circular geocentric 
orbits strongly perturbed by the sun and the moon. Programs based on 
them furthermore have shown that significant savings in computer time 
could be realized for the same accuracy, compared to the more common 
versions of the method of variation of parameters, with time as the in 


dependent variable. 
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APPENDIX A-2 


Auxiliary Developments in Singularity-Free Methods 


_ 2m 

A. 1 Mixed Fourier-Chebychev Series for p 
Since 

2m ,r . 2m 
r P 

we write, in the unperturbed motion 

- = 1 + (£*■ - l)cos E + i(~) c sin E 
a a a ah 

Since (-) 2 is a Fourier series in E, so will (~\) 2m using the 
a a 

recurrence 


r p r p p 


Finally, p Z is obtained as a Fourier-Chebychev series by multi- 

„ , . ,r.2m , ,a v2m 

plication of the senes for (-) and (~y 


A. 2 Mixed Fourier-Chebychev series for p cos 8 

Since, + + + 

r ' r i o 2 

p cos 0 = ~ a = (- ‘ 

r P a rr p a a r p 

r 

—2- has components which can be developed in Chebychev series 
a 

using the method of special values, with x as independent 
variable. - has components which are Fourier series in E, 


since 

£ = - + (r*- — ) cos E + -(~f)oSin E 
a a a a' a dE 

Therefore, p cos 0 will be a mixed Fourier-Chebychev series, 
a 2. 

because ( — ) is in turn a Chebychev series in x. 

r p 
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A. 3 Integration Algorithms 

[ 5 - 26 ] _ . 

Using the expressions given by Fox and Parker for the inte 

grals with respect to x of Chebychev's polynomials T r (x), and 
by a judicious use of integration by parts, the following algo- 
rithms, proved in Ref. [5-1], were obtained; 
dE 

(Note that ^ = constant = tt) 

K max 

a) Integration of ( R £^ SO(J,K)T R _^(x))sin jE 


Let the integral of the above expressions be; 


Km, 


Kmax 


( ! Sl(J,K)T K _ 1 (x))sin jE + ( R £ 1 Cl(J+l,K)T K _ 1 (x))cos jE 

K=1 


Then the algorithm follows, starting with the known SO's. 


* 0 


C1CJ + 1,K ) - - ] 

max J max 


S1(J,K -1) = - 

' ■’ ino v 


max 


2 (K - 1) 

ma - X -— ; Cl (J+1,K ) 

T /dE' max 

J V 


C1(J+1,K - 1) = ~ T S C (J,K - 1) 

max J max 


S1(J ,K - 2) = - 2 

’ max T db 

dx 


C1(J + 1, K - 1) 

x * tn o V 


max 


K 2 

Z2 = 2 — — S1(J,K - 2) 

max 

dE 

dx 


Zl= 


2 (K max 3) S1 (j k -2) 

_ . N 7 TV* OTT 


dE 


max 


dx 
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D0 1 K = K -3, 4,-1 
max 

C1(J + 1 ,K) = i(SO(J,K>- 22) 

Sl(J,K— 1) = -2 C1(J+1,K) + S1(J ,K+1) 

^ dx 

Z22 = Z2 

Z2 = Z1 

1 Zl = 2 K ~^ S 1 (J,K-1) + Z22 

dx 

C1(J+1, 3) = - j(S0(J,3) - Z2) 

S1(J , 2) ““ C^J+1,3) + SI (J,4) 

J dJ 

Z22 = Z2 

Z2 = Z1 

Z1 = ^ SI (J , 2) + \ Z22 

cx 

C 1 (J+1,2) =~i(S0(J,2) - Z2) 

S1(J,1) C1(J+1,2) + \ SI (J ,3) 

^ dx 

Z2 » Z1 

C1(J+1,1) = - 3 (S0(J,1) - z 2 ) 


END 
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K 

max 

b) Integration of { CO(J+l,K)T K _^(x) } cos jE 


Let the integral of the above expression be 

K K 

max max 

( r £ 1 S1 (J,K)T k _ 1 (x) ) sin jE + ( R £ 1 Cl(J+l,K)T K _ 1 (x))cos jE 


Algorithm is as follows, starting with the known SO T s 


C1(J+1, K ) = 0 
max 

S1(J, K ) = ~ C0(J+1, K ) 
max J max 


K .. 

C1(J+1, K - 1) = 2 ma ^" X S2(J,K ) 
v 5 max T dE max 

J dx 


K - 2 

C1(J+1, K - 2) » 2 ma ' X - r ~- S2 (J , K 

max T dE max 

J dx 


- 1 ) 


z2 = 2 ma * 
^ dE 

dx 


K - 2 

C1(J+1, K - 1) 

dE max 


K - 3 
max 


dx 


C1(J-KL, K - 2) 
max 


D0 


1 


K = K 

max 


3,4, -1 


S1(J,K) = - i(CO(J+l,K) - Z2) 

Cl (J+l, K-l) = 2 S1(J,K) + C1(J+1, K+l) 

^ dx 


Z22 = Z2 
Z2 = Z1 
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1 Z1 = 2 C1(J+1, K-l) + Z22 

dx 

S1(J,3) = 5 (C0(J+1, 3) -Z2) 

C1(J+1, 2) = ~~ S1(J,3) + C1(J+1, 4) 

Z22 = Z2 
Z2 = Z1 

Z1 = “ C1(J+1, 2) + i Z22 
dx 

S1(J , 2) = - J (C0(J+1, 2) -Z2) 

C1(J+1, 1) = S1(J ,2) + | C1(J+1, 3) 


Z2 = Z1 

1 

S1(J,1) = - (COCJ+1, 1) -Z2) 
c) Constant of integration 

The integrals being computed from E = 0 to E, we note that the 
perturbation should vanish at E = 0. Thus from the indefinite 
integral result, the value of the series at E = 0 £>r x = -1) 
should be subtracted i.e., in the final series (containing the 
constant of integration, J raax 

Cl(l,l) to be used = Cl(l,l) above - C1(J,K)T K _ 1 (-1) } 

the other coefficients remaining unchanged. 
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CHAPTER 6 

Orbital Programs 


6.1 Introduction 

In this chapter we briefly list, analyze and compare various 
programs, developed under this grant for the integration of 
perturbed orbits ^ 1 ^ , Other programs, such as SABAC, ECLIP, VOLER, 
etc. have been previously described. The program examined here 
differ from each other in the following aspects: 

a) the parameters, or osculating elements, being integrated 

b) the independent variable retained 

In order to make comparisons valud, it was decided to integrate 
by all methods examples which would serve as numerical standards in 
the analysis. These two reference solutions were obtained by NASA-' s 

numerical integration program ITEM, based on a modified Encke’s 

r 0-2] 

method and briefly described by Lowrey 
In all methods 

- identical integration techniques were used 

The same integration routine was used* a fourth order 
prediction-corrector method, of the Hamming type, with 
fourth-order Runge-Kutta-Gill starter* 
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- identical perturbation forces models 
These forces consist of 

1. "Third" body gravitational perturbation (by 
sun and moon) 

2. Gravitational perturbations due to the asphericity 
of the earth 

3. Atmospheric drag 

4. Solar radiation pressure 

The two last forces were set to zero on the examples treated. 
The position of the sun is interpolated from the American 
Ephemeris data for the sun's position at 0 h. E.T. every day. 
The moon is given by formulae having an error of + 0.75 min 
of arc at most, over 3 years from Jan 1, 1969. 

- identical reference standards 
6.2 Test 

The conditions for these are listed below, for two examples, one 
for close to zero eccentricity (e = 0.8 x 10 ), a = 0.58 x Earth- 

moon distance, orbital period = 12.05 days; one for high eccentricity 
(e = 0.936), a = 0.297 x Earth-moon distance, orbital period = 4.45 


days . 



INITIAL CONDITIONS TOR THE NUMERICAL EXAMPLES 


Example 1: Large Circular Orbit 

r = - 0.2932796 x 10 8 a 

O X 

- 2.0338597 x 10 8 a 

y 

8 '> 

+ 0.87225166 x 10 a meters 

z 

4 - 

(^) = 1.2405 x 10 3 a 

dt o x 

- 0.3358361 x 10 3 a 

y 

3 

- 0.36598403 x 10 a meters/second 

z 

T q = Feb. 18, 1971, 6.00 hours Ephemeris time 

From which 

r = 223235.8 K == 35 x radius of earth mean 
o m 

“ 0.58 x mean earth moon distance 

v = 1.336252 K /sec. 
o m 

a = 223234.0 K 

m 

e « 0.8018212 x 10 -5 
i = 28.50035° 

0 = 133.2179° 

u) = -54.99945° 
v = 180.0285° 

= 12.05 days 


Period 
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Example 2: Highly Eccentric Orbit 


r = -0.39275819084844 x 10 5 a 

O X 

-1.623139606007665 x 10 5 a y 

+0.8969904059411103 x 10 5 a K 

z m 

(4 1 ) = 0.1954377352706032 a 

dt o x 

-0.7854274357862668 a 

y 

+0.24190151060205798 a K /sec 

z m 

T q “ January 5, 1971, 18,5 hours Ephemeris time 

From which 

r = 189563,5 K 15 29*72 x radius of earth 
o m 

- 0,494 x mean earth-moon distance 


v 

o 

a 


0.8447536 K / sec 
m 

114151.4 K ~ 17.9 x radius of earth 
m 

- 0.297 x mean earth-moon distance 


e = 0.936227 
i - 33.40927° 
ft = 130.9163° 
co - -50.62353° 
v = 171.3767° 


Period 


4.45 days 
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The programs were designed to be "modular" in structure. Some 
modules, concerned with the computation of the perturbing forces, 
are identical with all the programs. Those dealing with integration 
are almost identical, except for the number of differential equations 
required and some print-outs. Such a modular arrangement makes it 
much easier to bring changes in some part of the computational scheme. 

A list of the programs to be discussed is given below 


NICE-T 

(for Numerical inte- 
gration of Circular 
& Elliptic Orbits) 

NICE-E 

NICE-EA 


BROUWER- E 
(Brouwer's method 
modified for circu- 
lar orbits) 

EOLA-T 


INDEPENDENT VARIABLE 
t 

E 

E; elements kept 
constant over (0,2ir) 
in E 

E 

t 


DEPENDENT VARIABLES 

-y 

r 0 * 

—>■ 

r O 5 

K i (i = 1,.. . ,6) 



Conventional oscu- 
lating elements 


6.2 Program NICE-T 

Independent variable: time 

Osculating parameters position vector r Q , and velocity vector r 0 , 
along the osculating orbit at time "T 0 " (fixed) 


Equations: six equations, in (5-2.60,61) of Chapter 5. 
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Tested: against NASA's ITEM, on the two examples quoted above. 

Thereafter used as a standard of comparison. 

Comparison with other programs: as most other programs had ”E" 

as independent variables at intervals of 2tt in E, 
data from these other programs consisting of time t; 
osculating parameters, elliptical osculating para- 
meters, radius and velocity vectors at t were punched 
out. At these times , the output of NICE-T was computed. 
The differences in the instantaneous values of the os- 
culating elements, r 0 and (— ) c , the elliptical para- 
meters and the values of the instantaneous radius and 
velocity vectors were computed and compared. The 
differences were then normalized by the maximum values 
of the quantities, as shown below. 

The integration spanned about 25 orbits, i.e. about 195 days in 
the high eccentricity case, and 300 days in the case of the large cir- 
cular orbit. 

Since their did not appear to exist any definite trend for these 
differences (except that the two first orbits had always much smaller 
differences than the remaining 23) , it was decided to "represent 15 them 
in the tables by their arithmetic mean A and standard deviation S„. 

A comparison of computer times is given in each case (CMU 360/67 
TSS) and will be commented upon for each program. 



RELATIVE ACCURACIES AND SPEEDS OF VARIOUS COMPUTER PROGRAMS 


Per unit errors in various quantities are defined as follows: 



£ i 180 

£2 — 0,j, 


360 

to — to T 
e w = 360 

V ' V T 
S 360 



where 

a = semi-major axis (meters)* 
e = eccentricity, 
i = inclination (degrees) 

ft = longitude of ascending nodes (degrees), 
u) = argument of perigee (degrees). 
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V = true anomaly (degrees). 

->■ 

r = instantaneous radius vector (meters), 
r = instantaneous velocity vector (meters/sec.). 

Quantities with subscript 'T* refer to values computed by program 
NICE T used as reference. 

A = arithmetic mean of 25 values (one at the end of each orbit for 25 
m 

orbits). 

Sp = standard deviation of the same 25 quantities. 

Example 1 - Large circular orbits (data as per example 1) 


Comparison of Computer Time/Orbit 
(Average of 25 Orbits) 



cpu Time/Orbit, Seconds 
16.9 

4.08 

2.8 

3.88 


BR0UWER E 


4.84 
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Example 2: Highly eccentric orbit (data as per example 2) 


Comparison of Computer Time/Orbit 
(Average of 25 Orbits) 


Program 

cpu Time/Orbit (Seconds) 

NICE T 

26.1 

NICE E 

5.48 

NICE EA 

3.07 

NICE EP 

5.21 

BROUWER E 

5.15 


E0LA - T _ 


9.95 



EXAMPLE 1 


COMPARISON OF ACCURACIES 



NICE 

E 

NICE EA 


NICE 

EP 

BROUWER 

E 


A 

in 

S D 

A 

m 

S D 

A 

m 

S D 

A 

m 

V 

e 

a 

0.13xl0' 6 

0.94x10" 12 

-0.18x10" 2 

0. 77xl0 -5 

0.81xl0~ 3 

0.59xl0" 6 

O.llxlO -5 

0.13xl0~ 10 

£ 

e 

0.65xl0 -6 

0.20x10" 8 

0.39xl0 -1 

0 . 19xl0 -2 

-.15x10 1 

.37xl0 -2 

0.19xl0" 4 

0.68xl0“ 7 

£ . 
1 

-0.27xl(f 6 

0.42xl0 -13 

-0.13xl0 -3 

O.llxlO" 6 

-0.40xl0 -3 

0.55xl0 -7 

-0 . 85xl0 -6 

-12 

0.62x10 

£ n 

-O.lOxlO -6 

-13 

0.12x10 

0, 57xl0~ 3 

O.llxlO" 6 

-0.17xl0" 3 

0.94xl0 -8 

-0. 70xl0 -7 

-12 

0.17x10 

£ 

U> 

-0.75xl0 -5 

0.45xl0 -10 

-0. 37xl0~ 2 

0. 26xl0 -2 

-O.lOxlO" 1 

0.79xl0 -3 

0.91xl0~ 5 

0.20xl0 -8 

e 

V 

0.94xl0" 5 

-10 

0.67x10 ±U 

-0. 74xl0" 2 

0.2 8x10" 2 

0. 20xl0~ 2 

0.74xl0 -3 

-0. 42xlO -4 

0. 26x1 0” 8 

e-> 

r 

-O.lOxlO -5 

0.39xl0 -11 

0.13xl0 _1 

0.34x10 4 

-0. 44xl0 -2 

0.16xl0 -4 

O.llxlO" 4 

0.74x10" 10 

o 

e-*~ 

r 

O.llxlO -5 

0.27xl0 -11 

-0. 14x10" 1 

0,53xl0" 4 

0.5xl0 -2 

0.18x10 4 

-O.llxlO" 4 

0.94xl0 -10 


i 


o 



EXAMPLE 2 

COMPARISON OF ACCURACIES 



NICE 

E 

NICE 

EA 

NICE 

EP 

BR0UWER E 

... ^ 

T 


A 

s 

A 

S„ 

A 

S^ 

A 

S,, 

A 

S n 


m 

D 

m 

D 

m 

D 

■ m 

D 

ra 

D 

£ 

a 

0.79xl0 -6 

0.32xl0 -10 

-0,46xl0 -2 

0.32xl0 -5 

0.12x10" 3 

0.34x10" 7 

0.56xl0" 5 

0.41xl0" 10 

-0. 19xl0 -4 

0.47xl0 -9 

£ 

e 

-0.74xl0" 6 

0.18xl0 -12 

-0.19xl0~ 3 

0.l4xl0 -7 

0. 28xl0~ 4 

0.72xl0" 9 

-o.eoxio" 5 

0. 17xl0~ 10 

-0.18xl0~ 4 

0.19xl0~ 9 

£ i 

0.85xl0“ 6 

-12 

0.36x10 

-0.13xl0 -3 

0.16x10" 7 

-0. 12xl0" 4 

0.18x10 9 

-0.13xl0 -5 

O.lOxlO" 10 

0. 29xlO -4 

0. 78xl0 -9 


0.34xl0 -6 

0.18xl0 -12 

- -0.90xl0 -4 

0,26x10 ^ 

-0.88xl0" 5 

0.51xl0" 9 

0.14xl0 -5 

O.lOxlO" 10 

-0. 20xl0~ 4 

0.46xl0~ 9 

£ 

CO 

- 

- 

-0.12xl0~ 2 

0.48xl0 -6 

-0. 66xl0~ 4 

0 . 87xl0" 9 

0.17xl0" 5 

0.25xl0 -11 

0.71xl0 -5 

0.15xl0 -9 

£ 

V 

0.45xl0" 4 

O.llxlO" 10 

0.85xl0 -2 

0.44x10 4 

-0. 79xl0 -3 

0. 6xl0~ 6 

0.76xl0“ 5 

O.lOxlO" 9 

O.llxlO" 3 

0. 14xl0~ 6 

£+ 

r 

0.46xl0 -4 

O.lOxlO -8 

0.56x10 ^ 

0.14xl0 -2 

0.74x10" 2 

0. 53xl0 -4 

0.57xl0" 4 

0.67x10 8 

0.64xl0~ 3 

0. 16xlO~ 4 

o 

e+ 

r 

-0. 14xl0 -3 

O.llxlO -7 

-0.36 

0.96xl0 -1 

0.24x10 ^ 

0 . 56x10" 3 

0. 22xl0~ 3 

0.98xl0 -7 


- 


ON 

I 



6-12 


6.3 NICE-E 

Independent variable: "E" of Chapter 5 (not eccentric anomaly) 
Osculating parameters: r Q > r 0 at E - 0 

Number of differential equations: 7 , given in 

of Chapter 5. 

Time required and accuracy: the method is about four times 

faster than that based on time as independent var- 
iable. At high eccentricities, a significant part 
of the gain could be that Keplers equation of time 
does not have to be inverted. However, the small 
time needed to do this inversion at small eccen- 
tricities fails to explain that the same gain still 
exists. Maybe the reason is to be found in the 
regularizing effect of using variable E, rather than 
time, which could be equivalently be seen as using a 
variable epoch time T 0 along the unperturbed orbit, 
in a way which presumably reduces computer time for 
a prescribed relative accuracy by keeping the motion 
in space (r 0 » r c , T c ) small. 


6.4 NICE-EA 

Independent variable: 


E modified, defined as -r= in 

/a 

x defined by ■— = 4— 4“ 

dx vp dr 


Osculating parameters: 

E 

Approximation mode:^— 


da 

dE 


r D , r 0 at E = 0 

<< 1 used in the definition of the 


independent variable; see Section 5. 
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Time required and accuracy: least of all methods* If the 

above approximation can be justified, either by the 
nature of the perturbing forces or, a posteriori by 
numerical experiments, this method could give a 20% 
saving in analytical integration (6 equations instead 
of 7 are used). The (relatively low) accuracy is 
of the same order as that obtained by a first-order 
Picard scheme. 

6.5 NICE-EP 

If the differential equations of motion are integrated in the 
iterations of a Picard scheme, no adequate analysis seems to 
exist on how long an interval of integration and how high an 
order of iteration one should take to minimize the computer 
time required for computing over a given interval with a pre- 
scribed error bound. 

Therefore, it was thought to be appropriate, to get an idea 
of what error is introduced only by retaining a first- order 
iteration scheme. To that effect, the equations of NICE-E 
were integrated over the typical (0, 2ir) interval for E while 
keeping the values of the elements constant over the interval. 

6.8 EOLA-NU 

This program integrates the classical elements a, e, i, m, £2 with 
respect to v, and obtains time t by integrating equation (2.3-4). 

No mean anomaly at epoch is used . The speed is comparable to that of 
E0LA-T . 
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CHAPTER 7 

General Conclusions 

At the conclusion of this work, we wish to briefly review the 
material developed under this grant and to make a few recommendations 
for future topics of study. 

New methods have been developed for the mission analysis and or- 
bital studies of satellites of the IMP-type, which describe trajectories 
strongly perturbed by the gravitational fields of the Sun and the Moon. 
These techniques cover a wide interval of the "accuracy vs. computer 
time" scale, ranging all the way from very fast methods of relatively 
low accuracy (as in SABAC) to high accuracy, more time-consuming schemes 
(as in the NICE programs) through a method based on non-numeric compu- 
tation, giving results of intermediate accuracy and time-consumption 
(as in VOLER) . All of these will be chosen at some point in the 
mission analysis :in the preliminary phase, in establishing large num- 
bers of possible launch windows; later on, in more detailed studies 
of better accuracy, and finally, in a few calculations by means of 
high accuracy programs suitable for low or high eccentricities, and 
which appear to save computer itme by a factor of 3 to 4, compared 
to conventional methods. 

Without going in detail into possible ways of implementing 
these suggestions, we think that various topics of investigation de- 
serve further study. One is the comparative analysis of existing 
methods, or the development of other techniques, which are most 
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suitable for multiple satellites, the dimensionality of the state 
vector (r c , r c ), say, going from 6 to 12, 18 etc* Another topic 
would be the development of literal theories based on regularized 
variables, since much insight is gained in the qualitative behavior 
of orbits even by means of "low-order 11 theories* Another area where 
more study appears to be needed is in minimizing computer time for a 
given error bound, or conversely, on a small computer of limited 
memory, in determining which method, for a given calculation time, 
assures the best accuracy in orbit determination* 



